How to Work With Your Own DNA Data

By the end of this you will have a local, indexed and reference-aligned copy of your own genotypes. It consists of a bgzipped VCF with a tabix index and a per-variant annotation table. The table holds gene, consequence, ClinVar assertion and gnomAD frequency. You will also have a QC report you can point at when a result looks strange, and a shell script that reproduces all of it from the original file. You need the original download from whoever sequenced or genotyped you, 200 GB of free disk if you have a CRAM or FASTQ (20 GB is plenty for array data), 16 GB of RAM, and a Unix shell. Everything below uses open tools: bcftools, samtools, plink2, bwa-mem2 or dragen-os, Ensembl VEP, and tabix. Install them with conda (mamba create -n dna -c bioconda -c conda-forge bcftools samtools plink2 bwa-mem2 ensembl-vep tabix) or use the official Docker images. This is measurement and interpretation. Nothing here is medical advice, and any variant you might act on needs a clinical geneticist and an orthogonal confirmatory test.

1. Identify which file you have

The phrase “DNA data” covers five very different things, and the first mistake people make is treating one like another.

Array export (23andMe, AncestryDNA, MyHeritage). A tab-separated text file. It is 15–25 MB uncompressed and holds around 600,000–700,000 rows. Columns are rsID, chromosome, position and genotype. This is a genotyping array: it reports the alleles at pre-chosen sites and says nothing about the other 3.1 billion bases. The site list is derived from SNP catalogs that began with efforts like the 1.42-million-SNP map published in 2001 ¹, and modern arrays still concentrate on common variation plus a curated set of known pathogenic sites.

VCF or gVCF from a sequencing provider. A single-sample VCF from a 30x whole genome is roughly 1–2 GB bgzipped and holds 4–5 million variants. A gVCF (with reference blocks) is 5–15 GB and lets you distinguish “matches reference” from “not called”.

CRAM or BAM. The aligned reads. A 30x CRAM against GRCh38 is 15–25 GB, and the equivalent BAM is 45–90 GB. This is the most useful file to own because every downstream call can be regenerated from it.

FASTQ. Raw reads, 2 x 50–60 GB gzipped for a 30x genome. Only worth keeping if you want to re-align against a different reference or a graph reference later.

Exome or panel data. Same formats, 1–2% of the genome. Check the BED file of targeted regions before concluding anything from an absent call.

Run this first, whatever you have:

# Array text
head -30 genome.txt
grep -c -v '^#' genome.txt

# VCF
bcftools view -h your.vcf.gz | grep -E '^##(reference|source|GATKCommandLine|DRAGEN)'
bcftools index -s your.vcf.gz | head          # variants per contig

# CRAM/BAM
samtools view -H your.cram | grep -E '^@RG|^@PG|^@SQ' | head -20

You are looking for three facts: the reference build, the caller, and the contig naming style (chr1 vs 1). Nearly every downstream failure traces back to one of these.

2. Pin down the reference build

Array files from consumer companies are almost always GRCh37/hg19 with 1-style contig names. Clinical and research sequencing since about 2018 is GRCh38 with chr1-style names. Mixing them produces coordinates that are silently wrong by tens to thousands of bases, which means the wrong gene, the wrong consequence, and an annotation table that looks plausible and is garbage.

Get a proper analysis set FASTA rather than a UCSC download:

# GRCh38, no-alt analysis set
wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa
samtools faidx GRCh38_full_analysis_set_plus_decoy_hla.fa

If your data is on GRCh37 and you want GRCh38, prefer re-alignment from reads over coordinate lifting. CrossMap/liftOver will drop or mis-place a few percent of sites, and the losses are concentrated in exactly the repetitive and structurally variable regions where you would most want accuracy. If you only have array text, lifting is your only option, so use CrossMap vcf with the UCSC chain file and keep the unmapped output file so you know what you lost.

CrossMap vcf hg19ToHg38.over.chain.gz array.b37.vcf.gz GRCh38.fa array.b38.vcf
# CrossMap writes array.b38.vcf.unmap - read it, do not delete it

3. Convert array text to a normalized VCF

Do not hand-parse the 23andMe file. Two properties will bite you. First, no-calls are written as -- or II/DD insertion-deletion codes, and a naive parser turns them into homozygous reference. Second, array genotypes are reported on the plus strand of the probe design, which is not necessarily relative to the reference allele. As a result, about half the time the “A/A” you see is the reference and half the time it is the alternate, and sometimes neither.

bcftools convert --tsv2vcf handles the reference-allele lookup for you:

# Strip comments, keep 4 columns, drop no-calls
grep -v '^#' genome.txt \
  | awk -F'\t' '$4 !~ /-|I|D/ {print $2"\t"$3"\t"$1"\t"$4}' \
  | sort -k1,1 -k2,2n > array.tsv

bcftools convert --tsv2vcf array.tsv \
  -f human_g1k_v37.fasta -s ME -Oz -o array.b37.vcf.gz
bcftools index -t array.b37.vcf.gz

Then normalize and check that the reference alleles agree with the FASTA:

bcftools norm -f human_g1k_v37.fasta -c ws -m -any \
  array.b37.vcf.gz -Oz -o array.norm.vcf.gz

-c ws warns about mismatches and sets the site to the FASTA allele. Count the warnings. More than a few thousand means your build guess is wrong. Rerun against the other build and compare.

For anything population-genetic (PCA, kinship, polygenic scores), convert to PLINK and work there instead:

plink2 --vcf array.norm.vcf.gz --make-pgen --out me
plink2 --pfile me --freq --missing --het --out me_qc

4. QC before you interpret anything

Every claim you make later rests on these numbers, so produce them once and keep them in the repo.

bcftools stats -F GRCh38.fa your.vcf.gz > stats.txt
grep -E '^SN|^TSTV' stats.txt

For a 30x single genome on GRCh38, expect roughly 4.0–5.0 million variants, 3.3–3.7 million SNVs plus 500–900k indels, a Ts/Tv ratio of 2.0–2.1 genome-wide, and 10,000–12,000 protein-altering variants. A Ts/Tv below 1.9 usually means a permissive filter or contamination. A variant count under 3.5 million usually means low coverage or aggressive filtering. Counts substantially above 5.5 million usually mean sample contamination or non-European ancestry relative to the reference (the reference bias is real, and rare-variant discovery in under-represented populations is systematically higher than the European baseline) ².

From a CRAM, check coverage and contamination directly:

samtools coverage your.cram | awk '$1 ~ /^chr[0-9XY]+$/'
mosdepth --by 1000 --fast-mode me your.cram     # then read me.mosdepth.summary.txt
verifybamid2 --SVDPrefix 1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38.fa --BamFile your.cram

FREEMIX above 0.02 means measurable contamination and you should not trust heterozygous calls at low depth. Also confirm the sex chromosomes match expectation: chrX at half autosomal depth with chrY present, or chrX at full depth with chrY near zero. A mismatch between the file and your expectation means a sample swap, which happens more often than labs like to admit, and it invalidates everything downstream.

5. Annotate with VEP

Annotation is where you go from coordinates to meaning. Use Ensembl VEP with a local cache so you are not rate-limited and your results are reproducible.

vep -i your.vcf.gz -o annotated.vcf.gz \
  --cache --dir_cache $HOME/.vep --offline \
  --assembly GRCh38 --fasta GRCh38.fa \
  --vcf --compress_output bgzip \
  --everything \
  --pick_allele_gene --pick_order canonical,rank,length \
  --custom file=clinvar_20250601.vcf.gz,short_name=ClinVar,format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN \
  --plugin CADD,snv=whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz \
  --fork 8 --stats_file vep.html

Notes on the flags that matter. The --pick_allele_gene flag gives you one consequence per gene per allele instead of twenty transcript rows, which makes the output tractable. The --offline flag guarantees the annotation version is whatever is in your cache, so record the cache release number. Setting type=exact on the ClinVar custom file makes the match position-and-allele specific rather than positional. That prevents you from inheriting a pathogenic label from a different allele at the same coordinate. The CLNREVSTAT field is as important as CLNSIG. A “Pathogenic” call with review status no_assertion_criteria_provided from a single submitter is weak evidence, and a large share of ClinVar entries carry conflicting or low-confidence assertions.

Flatten to a table you can grep:

bcftools +split-vep annotated.vcf.gz -d -f \
  '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Consequence\t%gnomADg_AF\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT\t%CADD_PHRED\t[%GT]\n' \
  -i 'GT!="0/0" && GT!="./."' > annotated.tsv

A useful triage filter: consequence in the high-impact set, gnomAD genome allele frequency below 0.001, and either a ClinVar assertion with at least one-star review status or a CADD PHRED above 20. On a typical genome that leaves a few hundred rows. Most are benign despite the filter, because rarity and predicted impact together are still weak predictors of consequence in an individual.

Expect the bulk of your variants to be noncoding, and expect interpretation there to be thin. The ENCODE pilot showed that functional elements blanket the genome well beyond protein-coding exons, including large numbers of transcribed and regulatory regions with no obvious phenotype attached ³. Twenty years later, a noncoding variant with a strong regulatory annotation and no disease association is still uninterpretable for a single individual.

6. Answer a specific question

Targeted queries are where personal genome data is genuinely useful, because you are checking a small set of well-characterized sites rather than fishing.

# One gene, GRCh38 coordinates from Ensembl
bcftools view -r chr13:32315474-32400266 annotated.vcf.gz | bcftools +split-vep -f '...'

# A specific rsID, if your VCF is annotated with dbSNP IDs
bcftools view -i 'ID="rs1801133"' annotated.vcf.gz

Pharmacogenes are a special case. Most clinically meaningful pharmacogenetic results are star alleles, which are haplotypes defined by combinations of variants plus gene copy number, and CYP2D6 in particular has deletions, duplications, and hybrid genes with CYP2D7 that a per-variant VCF cannot represent. Use PharmCAT rather than reading individual positions:

java -jar pharmcat-pipeline.jar -vcf your.vcf.gz -o pharmcat_out

PharmCAT emits diplotype calls and the CPIC guideline text. Read it as a measurement. The output is not an instruction. Drug selection and dosing are clinical decisions and require a physician or pharmacist who can see your full picture.

7. Know what the file cannot tell you

Array data reports only the designed sites ¹. If a variant is not on the chip, your file has no information about it, and the absence of a row is not the absence of a variant. Imputation against a reference panel recovers common variation well and rare variation poorly, and panel coverage is uneven across ancestries, so imputed rare genotypes in an under-represented population are unreliable ².

Short-read whole genome sequencing, even at 30x, is blind or weak in several places. Homopolymers and short tandem repeats are miscalled, so expansion disorders need a targeted assay. Segmental duplications (SMN1/SMN2, PMS2, the CYP2D6 region) give ambiguous alignments. Copy number changes below a few kilobases are detectable only with dedicated callers like manta or cnvkit, and those callers have high false-positive rates on single samples. HLA typing needs HLA-LA or a similar graph-based tool. The standard pipeline is not enough for it. Low-level somatic mosaicism below roughly 10% allele fraction is indistinguishable from sequencing error at 30x, because the per-base error floor of reversible-terminator chemistry sits in that same range ⁴. And nothing in a germline genome tells you about methylation, expression, or current physiology.

8. Store it, and decide what you share

Keep three copies: working disk, a local archive, and an encrypted offsite copy. Archive the CRAM and the reference FASTA together, since a CRAM without its exact reference is unreadable.

# Encrypt before upload
age -r age1... -o dna.cram.age your.cram
# or
gpg --encrypt --recipient you@example.com your.cram
sha256sum your.cram > your.cram.sha256

On sharing: uploading a raw array file to a free analysis site or a genealogy matching database is an irreversible disclosure. Your genotypes identify your relatives as well as you, which is the mechanism behind forensic genetic genealogy, where investigators build family trees from public match databases to identify people who never uploaded anything. The United States operates CODIS, a national forensic database of STR profiles from arrestees, convicted offenders, and crime scenes. It is not a database of everyone, and it does not hold sequence data. But no one needs your sample to find you if a second cousin has uploaded theirs. Consumer genetics companies are also corporate assets that get sold, breached, and subpoenaed, and their privacy policies are revisable.

In the US, GINA restricts the use of genetic information in health insurance and employment. It does not cover life, disability, or long-term care insurance. The gap is the practical reason to think carefully before putting identified genetic data anywhere you do not control.

Common problems

Coordinates look right but the gene is wrong. Build mismatch. Check ##reference in the header and spot-check a known SNP: rs1801133 is at chr1:11856378 on GRCh37 and chr1:11796321 on GRCh38.

bcftools reports “the sequence not present” or zero overlaps. Contig naming mismatch. Fix with bcftools annotate --rename-chrs chr_map.txt where the map is two columns, old and new.

Reference-allele mismatches after bcftools norm -c ws. For array data, the genotype was on the opposite strand or the site is a known strand-ambiguous A/T or C/G SNP. Drop strand-ambiguous sites from any analysis that depends on allele identity rather than accepting a guess.

Your VCF says homozygous reference at a site you know you carry. Check whether it is a real call or an untested position. A single-sample VCF omits reference-matching positions entirely, so absence means nothing. Use the gVCF or compute depth from the CRAM at that exact coordinate with samtools depth -a -r chr1:11796321-11796321.

VEP output has thirty rows per variant. You omitted --pick_allele_gene, or you are reading the per-transcript output as if it were per-variant.

A polygenic score comes out at an extreme percentile. Scores are calibrated in the ancestry they were trained in and shift substantially across populations. Check the reference distribution the score was derived from before believing the percentile.

A ClinVar “Pathogenic” hit in a gene you have never heard of. Check CLNREVSTAT. Check the gnomAD frequency (a truly pathogenic dominant variant will not be at 1% in gnomAD). Check read support in the CRAM with samtools tview, and then take it to a clinical geneticist for confirmatory testing. Single-sample short-read calls at clinically relevant sites carry a false-positive rate high enough that no one should act on one unconfirmed.

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

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

4. David Bentley, Shankar Balasubramanian, Harold Swerdlow, et al. Accurate whole human genome sequencing using reversible terminator chemistry. Nature, 2008. https://doi.org/10.1038/nature07517 ↩