How to Download Your DNA Data and Turn It Into Something You Can Analyze

Oak
genomics vcf plink imputation raw-dna
A lab instrument pulls a translucent DNA ribbon with sparse glowing beads through a glass column while blue mist fills the gaps with fainter beads.

By the end of this guide you will have your consumer genotype file downloaded and checksummed, converted into a sorted, bgzipped, reference-correct VCF on GRCh38 with an index, optionally imputed against a modern reference panel, and annotated so you can query it with bcftools or hand it to an agent. You need a shell, roughly 20 GB of free disk if you plan to impute, and four tools: plink2, bcftools (1.17 or newer), samtools, and either CrossMap or Picard for the liftover. Everything here runs fine on a laptop except imputation, which we will offload to a server. The whole path, download to annotated VCF, takes about an hour of attention and a few hours of waiting.

One framing point before the steps. A consumer DNA test is a genotyping array: it interrogates a fixed set of pre-chosen sites, typically 600,000 to 700,000 of them on Illumina’s Global Screening Array, and reports which of two expected alleles you carry at each. The human genome is about 3.1 billion bases 1, so the array reads roughly 0.02% of your sequence. That is enough for ancestry, relative matching, and polygenic scores built on common variation. It is not sequencing, and the difference matters for almost every clinical question, which is the last section of this guide.

1. Download the raw file from the vendor

Each vendor buries the export in account settings, and each hands you a slightly different format.

On Ancestry, go to your DNA test’s settings page (DNA → Your DNA Results Summary → Settings), find “Download Raw DNA Data,” confirm by email, and wait for a link. You get a ZIP containing a single tab-separated AncestryDNA.txt. On 23andMe, use Profile → Settings → scroll to “23andMe Data” → Download raw data; you receive genome_<name>_Full_<date>.zip containing a tab-separated file with a long comment header. MyHeritage puts it under Manage DNA kits → the three-dot menu → Download raw data, and emails a ZIP with a CSV. FamilyTreeDNA offers an autosomal “Build 37 Raw Data Concatenated” CSV under Data Download.

To answer the question that gets asked most often: there is no legitimate way to get your own genome sequenced for free. The free options are (a) exporting data you already paid for, which is what this guide is about, (b) research cohorts that return some results in exchange for enrolling and contributing your data, and (c) analyzing public data such as 1000 Genomes or Personal Genome Project releases, which are other people’s genomes and useful only for building and testing pipelines. If someone offers a free “DNA upload analysis,” you are paying with the file.

Checksum the archive the moment it lands, because truncated downloads are the single most common failure here and they fail silently at the parsing stage:

unzip -p AncestryDNA.zip > AncestryDNA.txt
shasum -a 256 AncestryDNA.txt | tee AncestryDNA.txt.sha256
wc -l AncestryDNA.txt

Expect 600,000 to 720,000 data lines for a recent Ancestry or 23andMe v5 kit, around 960,000 for 23andMe v3, and about 700,000 for MyHeritage. Anything under 100,000 lines means the file is truncated or you downloaded a filtered subset.

2. Read the header and establish the genome build

Before any conversion, look at what you have. This is 30 seconds that prevents an entire bad analysis.

head -25 AncestryDNA.txt

Ancestry’s header states the array version and, importantly, that positions are on GRCh37 (build 37, hg19). The data columns are rsid, chromosome, position, allele1, allele2, and chromosomes are coded numerically with 23 = X, 24 = Y, 25 = pseudoautosomal region of X, 26 = mitochondrial. No-calls are 0. 23andMe’s file also uses build 37, gives a single concatenated genotype column (AA, AG, -- for a no-call), and codes the sex chromosomes as X, Y, MT. Both vendors still ship build 37 as of this writing, so assume 37 unless the header says otherwise.

Check your overall call rate, since a poor-quality sample shows up here first:

awk 'BEGIN{FS="\t"} !/^#/ && NR>1 {t++; if ($4=="0"||$5=="0") n++} END{printf "sites=%d nocall=%d rate=%.4f\n", t, n, 1-n/t}' AncestryDNA.txt

A healthy kit lands above 0.98. Below about 0.95, the vendor’s own QC usually would have flagged it, and imputation quality degrades noticeably.

3. Normalize to a single input format and convert to VCF

plink2 reads the 23andMe layout directly with --23file, so the cleanest route for Ancestry and MyHeritage files is a small awk step that produces that layout. For Ancestry:

awk 'BEGIN{FS=OFS="\t"}
  /^#/ {next}
  $1=="rsid" {next}
  {
    chr=$2
    if (chr==23) chr="X"; else if (chr==24) chr="Y";
    else if (chr==25) chr="XY"; else if (chr==26) chr="MT";
    g=$4 $5
    gsub(/0/,"-",g)
    print $1, chr, $3, g
  }' AncestryDNA.txt > input23.txt

Then convert. Two passes, because --ref-from-fa needs sorted variants:

plink2 --23file input23.txt MYID --snps-only just-acgt \
       --sort-vars --make-pgen --out step1

plink2 --pfile step1 --fa GRCh37.primary.fa --ref-from-fa force \
       --output-chr chrM --export vcf bgz id-paste=iid --out mydna_b37

tabix -p vcf mydna_b37.vcf.gz

--snps-only just-acgt drops indel-like and non-ACGT probes that arrays report inconsistently. --ref-from-fa force is the step everyone skips and should not: array files do not tell you which allele is the reference, and plink2 will otherwise pick the major allele, producing a VCF where REF and ALT are wrong at a large fraction of sites. Downstream annotation then silently returns garbage. You need a build 37 FASTA whose contig names match your --output-chr choice; Ensembl’s primary assembly FASTA is a reasonable source, and Ensembl is also where rsID-to-coordinate mapping and transcript annotation ultimately come from for most tools you will use 2.

Verify the reference match rate immediately:

bcftools norm -c w -f GRCh37.primary.fa mydna_b37.vcf.gz -Ou -o /dev/null

-c w warns instead of failing. You want zero or near-zero REF mismatches. Thousands of mismatches means your FASTA is the wrong build, which brings us to the next step.

4. Lift over to GRCh38

Almost every current reference resource (gnomAD v4, ClinVar’s primary VCF, the TOPMed imputation panel, VEP’s default cache) is on GRCh38. Working on build 37 means constantly translating coordinates, and coordinate translation is where mistakes enter. Lift once, early.

CrossMap vcf GRCh37_to_GRCh38.chain.gz mydna_b37.vcf.gz \
         GRCh38.primary.fa mydna_b38.vcf

bcftools sort mydna_b38.vcf -Oz -o mydna_b38.sorted.vcf.gz
tabix -p vcf mydna_b38.sorted.vcf.gz
bcftools norm -c w -f GRCh38.primary.fa mydna_b38.sorted.vcf.gz -Ou -o /dev/null

CrossMap will swap alleles where the reference changed between builds and write the rest to an .unmap file. Expect to lose a few thousand sites out of 650,000, concentrated in regions that were restructured between builds. Picard’s LiftoverVcf with RECOVER_SWAPPED_REF_ALT=true and a REJECT output is an equally good choice and gives a more detailed rejection log. Read that log rather than deleting it: a cluster of rejections in a region you care about is information.

5. Impute, if you want more than the probe set

Imputation uses the correlation structure of haplotypes in a reference panel to infer genotypes at sites your array never measured. Done well, it takes you from 650,000 genotyped sites to roughly 25 to 40 million imputed sites with usable confidence. It is the highest-value single operation you can perform on an array file, and it is also where people most often overinterpret the output.

Use the TOPMed Imputation Server rather than running Beagle locally against 1000 Genomes. The TOPMed panel is much larger and imputes low-frequency variation far better, and the server handles phasing (Eagle) and imputation (minimac4) for you. Prepare per-chromosome, bgzipped, sorted, GRCh38 VCFs with chr-prefixed contig names:

for c in $(seq 1 22); do
  bcftools view -r chr${c} mydna_b38.sorted.vcf.gz -Oz -o chr${c}.vcf.gz
  tabix -p vcf chr${c}.vcf.gz
done

Submit with build GRCh38, phasing Eagle, and population “all” for QC. The server’s QC report tells you how many sites matched the panel, how many were strand-flipped, and how many were excluded for allele mismatch. A normal array gives over 95% panel overlap. If your overlap is below 80%, your REF/ALT assignment in step 3 was wrong.

Then filter the results. Each imputed site carries an R2 INFO field, minimac’s estimate of the squared correlation between imputed and true genotype dosage:

bcftools concat -Oz -o imputed.vcf.gz chr*.dose.vcf.gz
bcftools +fill-tags imputed.vcf.gz -Oz -o imputed.tagged.vcf.gz -- -t MAF
bcftools view -i 'INFO/R2>0.8 && INFO/MAF>0.01' imputed.tagged.vcf.gz \
         -Oz -o imputed.common.vcf.gz

We use R2 > 0.8 for common variants in polygenic score work and would not trust a single imputed site below about R2 > 0.9 for anything you intend to look at individually. The structural limit is that imputation cannot recover variation absent from the panel: a private or family-specific variant has no haplotype neighbors to borrow from and will never appear, no matter how good the panel is. This is the reason imputed array data does not substitute for sequencing when the question concerns rare variants.

6. Annotate and query

Now the file becomes useful. Annotate consequences with VEP and cross-reference population frequencies:

vep --input_file mydna_b38.sorted.vcf.gz --output_file annotated.vcf.gz \
    --vcf --compress_output bgzip --cache --assembly GRCh38 \
    --everything --check_existing --pick_allele_gene --fork 4

VEP’s transcript models and cross-references come from Ensembl, so your consequence calls are only as current as the cache release you install, and a variant’s “canonical transcript” can change between releases 2. Pin the release version in whatever you build on top of this.

For lookups, bcftools is faster than loading anything into a dataframe:

bcftools query() { bcftools view -H -r "$1" annotated.vcf.gz; }
bcftools view -H -i 'ID=="rs1801133"' annotated.vcf.gz | cut -f1-5,10

Two interpretation constraints belong here rather than at the end. First, array genotypes at individual clinically relevant sites carry a meaningful false-positive rate, because a single probe with a nearby polymorphism can miscall consistently, and consumer arrays are not run as diagnostic tests. A variant that appears in your raw file and looks alarming in ClinVar needs orthogonal confirmation by sequencing in an accredited clinical laboratory before it means anything, and the interpretation of a confirmed result belongs with a clinical geneticist or genetic counselor. Second, absence of a variant in an array file is close to meaningless: the array did not look at the vast majority of sites in any gene.

7. Know what the array file cannot contain

Understanding the boundaries of the file is part of working with it well. Your array VCF has no information about indels beyond a handful of probed sites, no copy-number or structural variants, no short tandem repeat lengths, and essentially no coverage of the genes and regions that array designers did not target. It has no RNA-level information at all, so nothing about which transcripts you express or at what level; that requires a separate assay, and the conceptual link between expressed sequence and genome annotation goes back to the earliest cDNA tag work 3.

Mitochondrial DNA deserves a specific warning. Arrays probe a few hundred mitochondrial sites at most, which is enough for haplogroup assignment and nothing else. Even with whole-genome sequencing, mitochondrial variant calling is complicated by nuclear-embedded mitochondrial DNA segments, which were found to be widespread across 66,083 sequenced human genomes and can masquerade as low-level mitochondrial heteroplasmy if reads are not carefully assigned 4. If you see a heteroplasmy claim from any pipeline, check how it handled NUMTs before believing it.

Common problems

Reference allele mismatches after liftover usually mean the wrong FASTA, not a broken chain file. Confirm with bcftools norm -c w against both builds and see which one matches; if neither does, your original conversion never applied --ref-from-fa.

Strand ambiguity at A/T and C/G sites is the classic array trap. These variants cannot be strand-resolved by allele identity alone, so if you merge your file with any external dataset you will introduce errors at roughly 15% of sites unless you resolve strand by allele frequency. The practical answer is to exclude ambiguous sites from any merge or to let the imputation server’s QC handle it, since the TOPMed pipeline checks strand against the panel by frequency and reports flips explicitly.

Duplicate and multi-allelic probe entries cause plink2 and bcftools to disagree about variant counts. Deduplicate by position and allele before conversion with plink2 --rm-dup force-first, and split multi-allelics with bcftools norm -m-any before annotation.

Heterozygous calls on chromosome Y or the non-pseudoautosomal X in a male sample indicate probe-level noise, not biology. Use plink2 --split-par b38 to separate the pseudoautosomal regions correctly and treat residual heterozygosity there as a quality flag on that region.

No-call runs clustered in one chromosome arm point to a sample or scan problem rather than a deletion. Compare the no-call density per 10 Mb window across the genome before interpreting any apparent gap:

bcftools query -f '%CHROM\t%POS\t[%GT]\n' mydna_b38.sorted.vcf.gz \
  | awk '$3=="./."{print $1"\t"int($2/1e7)}' | sort | uniq -c | sort -rn | head

Finally, if your goal in downloading the file was to answer a question about a specific gene, a specific rare variant, or anything you would take to a physician, an array export will not get you there, and neither will imputing it. That requires sequencing at depth, and the interpretation requires a clinician. The array file is excellent for ancestry, relative matching, polygenic scores over common variation, and learning a genomics toolchain on your own data, which is a genuinely good reason to do all of the above.

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. Tabitha M Powledge. Human genome project completed. Genome Biology, 2003. https://doi.org/10.1186/gb-spotlight-20030415-01 ↩

  2. T. Hubbard. The Ensembl genome database project. Nucleic Acids Research, 2002. https://doi.org/10.1093/nar/30.1.38 ↩ ↩2

  3. Mark D. Adams, Jenny M. Kelley, Jeannine D. Gocayne, et al. Complementary DNA Sequencing: Expressed Sequence Tags and Human Genome Project. Science, 1991. https://doi.org/10.1126/science.2047873 ↩

  4. Wei Wei, Katherine R. Schon, Greg Elgar, et al. Nuclear-embedded mitochondrial DNA sequences in 66,083 human genomes. Nature, 2022. https://doi.org/10.1038/s41586-022-05288-7 ↩