How to Analyze Your Raw DNA Data

This walkthrough is for anyone who has a consumer DNA file sitting on a hard drive and wants to do something rigorous with it. The goal is to take you from a vendor text file to a defensible set of results. By the end you will have a normalized, build-consistent VCF of your own genotypes, annotated against ClinVar and gnomAD. You will also have a short list of variants worth a conversation with a clinical geneticist, an ancestry and haplogroup call you can defend, and a clear sense of which questions your file cannot answer.

The prerequisites are modest. You need the raw export from a consumer array, which is a .txt or .zip from 23andMe, AncestryDNA or MyHeritage. FASTQ or CRAM files from a sequencing provider also work. On the hardware side you need a Linux machine with 16 GB RAM and 50 GB free for array work. If you are processing a 30x genome, plan on 32 GB RAM and roughly 500 GB of disk. You also need bcftools, samtools, plink2, Ensembl VEP, and Docker installed.

One framing note before we start. Everything below is measurement and interpretation. Nothing here is a diagnosis. Anything that touches disease risk or medication needs a clinician, specifically a board-certified medical geneticist or genetic counselor with access to a CLIA/CAP-certified confirmatory lab.

1. Get the file and understand what it contains

The first task is retrieving your raw data and forming an accurate mental model of what the vendor measured. This matters because most downstream mistakes trace back to a misunderstanding at this stage.

In 23andMe, the path is Account → Browse Raw Data → Download → “Request Download”. In AncestryDNA, it is Settings → Download Raw DNA Data, confirmed by email. Both vendors deliver a gzipped tab-separated file of roughly 15 MB compressed.

The 23andMe format consists of four columns following a block of # comments:

# rsid  chromosome  position  genotype
rs4477212   1   82154   AA
rs3094315   1   752566  AG
i713426     1   787173  --

Several properties of this file shape everything that follows. They are worth internalizing before you run a single command.

• This is a genotyping array rather than sequencing. The v5 chip interrogates roughly 630,000 sites, which is about 0.02% of the 3.1 Gb genome. Those sites were chosen for ancestry informativeness and common-variant coverage rather than clinical completeness.
• Positions are reported on GRCh37 (hg19) by every major consumer vendor to date. Nearly every current annotation resource is built on GRCh38, so you will be lifting over between the two.
• A genotype of -- means the probe failed to produce a call. In AncestryDNA files, I and D denote an insertion or deletion with no allele sequence attached. That makes those rows useless to downstream tools.
• Genotypes for some probes are reported on the chip’s design strand. A/T and C/G sites are therefore strand-ambiguous, and they are a common source of silent errors.
• Identifiers beginning with i are 23andMe-internal probes with no public rsID mapping.

2. Convert to VCF against a reference

The vendor text file is not a format any serious tool accepts. The next step is converting it into a VCF, the standard tabular format for describing variants relative to a reference genome. The conversion also gives us a chance to correct the file against that reference rather than trusting it.

bcftools convert --tsv2vcf handles the 23andMe layout directly. It also sets the REF allele from the reference FASTA rather than taking the vendor’s word for it, which matters more than it sounds:

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz
gunzip human_g1k_v37.fasta.gz && samtools faidx human_g1k_v37.fasta

zcat genome_raw.txt.gz | grep -v '^#' | awk '$4!="--"' > clean.tsv

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

Sites whose alleles do not match the reference are dropped with a warning. Read that warning count rather than scrolling past it. If more than about 1% of sites fail, your chromosome naming or your build assumption is wrong. The rest of the analysis will inherit the error.

With a valid GRCh37 VCF in hand, lift it to GRCh38. Then normalize the representation of each variant so that annotation tools match it correctly:

picard LiftoverVcf I=me.b37.vcf.gz O=me.b38.vcf.gz \
  CHAIN=hg19ToHg38.over.chain.gz R=GRCh38_full_analysis_set.fa \
  REJECT=rejected.vcf

bcftools norm -f GRCh38_full_analysis_set.fa -m -any -c w \
  -Oz -o me.norm.vcf.gz me.b38.vcf.gz
bcftools index me.norm.vcf.gz

Expect a few thousand rejected sites. Most of them sit in regions that were restructured between the two builds. Look through rejected.vcf for anything you care about rather than assuming the whole file is noise.

3. Quality-check before you interpret anything

Before any variant means anything, you need evidence that the file is intact and that it belongs to the person you think it does. The three checks below take about two minutes. They catch file mix-ups, sample swaps, and vendor chip changes.

plink2 --vcf me.norm.vcf.gz --make-bed --out me
plink2 --bfile me --missing --het --check-sex --out qc

Read the three outputs as follows.

• Call rate. qc.smiss should show missingness under 2%. Above 5%, the array run was marginal. Heterozygous calls are the first thing to become unreliable.
• Heterozygosity. The F statistic in qc.het should sit near 0 for an outbred individual. Strongly negative F means excess heterozygosity and usually indicates that two samples were mixed. Strongly positive F points to either consanguinity or allele dropout.
• Sex check. X heterozygosity should be near 0 for XY and around 0.3 for XX. A mismatch with what you expect means either the file is not yours or there is a sex-chromosome aneuploidy. The latter is a clinical finding rather than a bug to route around.

4. Annotate with VEP, ClinVar, and gnomAD

With a clean callset, you can attach biological meaning to each variant. Annotation tells you which gene and transcript a variant falls in, what ClinVar says about it, and how common it is in large population databases such as gnomAD.

Use the offline VEP cache for this. Online annotation services that ask you to upload your file are taking possession of a permanent identifier for you and your relatives. We return to that point in step 9.

docker run -v $HOME/vep:/data ensemblorg/ensembl-vep \
  vep --offline --cache --dir_cache /data --assembly GRCh38 \
      --fasta /data/GRCh38.fa --fork 8 --vcf --compress_output bgzip \
      -i /data/me.norm.vcf.gz -o /data/me.vep.vcf.gz \
      --everything \
      --custom /data/clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
      --custom /data/gnomad.genomes.v4.1.sites.chr%CHR%.vcf.bgz,gnomADg,vcf,exact,0,AF,AF_nfe,AF_afr,AF_eas

Now pull out anything ClinVar flags, carrying the review status along with the label. The review status tells you how much evidence stands behind a classification, and it matters more than the label itself:

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

Filter that list hard. Drop anything whose CLNREVSTAT is no_assertion_criteria_provided or criteria_provided,_single_submitter. Also drop anything with a gnomAD allele frequency above roughly 0.1% if the associated condition is a rare dominant disorder. ExAC demonstrated the scale of this problem. A large fraction of variants previously reported as disease-causing turn out to be present in population reference cohorts at frequencies far too high to explain a rare Mendelian phenotype ¹. Ancestry matters here as well, because allele frequencies estimated from a panel that does not include your population will mislead you. Deep sequencing of population-specific cohorts has repeatedly produced large numbers of rare variants absent from earlier global panels ².

5. Know why most “pathogenic” array hits are false

This is the single most important step in the whole exercise, and it is conceptual rather than computational. Understanding it will save you from the most common and most distressing error in consumer genomics.

An array probe is a hybridization assay with a fixed error rate per site. For a common variant present at 30% frequency, a 99.9% accurate call is perfectly serviceable. For a variant present in 1 in 50,000 people, the prior probability is so low that even a small per-site error rate makes most positive calls wrong. The arithmetic is unforgiving. If the true carrier rate is 2×10⁻⁵ and the false-positive rate is 10⁻³, the vast majority of flagged carriers are artifacts. Rare pathogenic variants are precisely the class where consumer arrays perform worst, and they are precisely what the free upload sites will surface for you.

The practical rule follows directly. Treat every rare variant call from an array as a hypothesis with maybe a 1-in-10 chance of being real. Confirm it by targeted Sanger or clinical-grade sequencing through a physician before it changes anything. Do not act on clinvar_hits.tsv. Take it to a genetic counselor.

The second limitation is coverage, and it cuts in the opposite direction. The array measures only the probes it carries, so an absent BRCA1 finding tells you nothing. The chip covers a handful of BRCA1 positions out of thousands of known pathogenic variants, and it covers no copy-number changes at all.

6. Ancestry and haplogroups, done properly

Ancestry inference is where array data is genuinely strong. It depends on common variants measured across many sites, which is exactly what the chip was designed to do.

Start with the uniparental lineages, which trace the direct paternal and maternal lines:

# Y haplogroup from the array's Y calls
yhaplo -i genome_raw.txt --format 23andMe -o yhaplo_out

For the mitochondrial haplogroup, extract the MT rows and submit them to HaploGrep’s local jar. Array coverage of the mitochondrial genome is partial, a few thousand probes at most. Expect a terminal clade one or two levels shallower than full mitochondrial sequencing would give you.

Autosomal ancestry requires merging your genotypes with a public reference panel. You then run ADMIXTURE or a principal components analysis:

plink2 --bfile me --bmerge ref_panel --geno 0.02 --maf 0.05 \
       --indep-pairwise 200 50 0.2 --make-bed --out merged
admixture --cv merged.bed 8 -j8

Choose a panel with real population depth. The Simons Genome Diversity Project provides 300 high-coverage genomes from 142 populations. It makes a far better backbone than the small, Eurocentric panels bundled with hobbyist tools ³. If deep-time components interest you, bear in mind that the ancestral composition of some groups is genuinely admixed at the source. The 24,000-year-old Mal’ta genome showed that Native American ancestry derives from a mixture of an ancient north Eurasian population and East Asians. That is why single-source assignment models produce incoherent output for some individuals ⁴.

Read the resulting percentages as distances to reference populations under an explicit model, rather than as facts about your ancestors.

7. Pharmacogenomics: what the array can and cannot call

Pharmacogenomics is the study of how genotype affects drug response. It is an instructive case of array data failing in a specific, predictable way. Running PharmCAT is a good way to see the structure of the problem:

java -jar pharmcat-pipeline.jar me.norm.vcf.gz -o pgx/

PharmCAT will report large numbers of uncalled star alleles from array input, and it is correct to do so. Diplotype calling requires that all defining positions for a haplotype be observed. A single missing position collapses the result to “no call” or to a default *1, which silently overstates normal enzyme function. CYP2D6 is worse still, because its major functional variation is structural. That includes whole-gene deletions, duplications, and CYP2D6/CYP2D7 hybrids that an array does not measure at all.

Read the output as a map of what is missing rather than as a report of what is true. Anything that touches a medication decision goes through a clinician and a clinical pharmacogenomics lab, not through your terminal.

8. If you have whole-genome sequencing, build the callset yourself

If your provider gave you FASTQ or CRAM files rather than a finished VCF, you control the entire pipeline from alignment through variant calling. That also means you can measure its performance directly. This is the version we would run:

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:L1\tSM:ME\tPL:ILLUMINA\tLB:lib1' \
  GRCh38_full_analysis_set.fa R1.fastq.gz R2.fastq.gz \
| samtools sort -@ 8 -m 4G -o me.bam -
samtools index me.bam

docker run -v "$PWD":/in google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS --ref=/in/GRCh38_full_analysis_set.fa \
  --reads=/in/me.bam --output_vcf=/in/me.dv.vcf.gz \
  --output_gvcf=/in/me.dv.g.vcf.gz --num_shards=32

Having produced a callset, measure it before you believe it:

mosdepth --by 1000 --fast-mode me.dist me.bam
bcftools stats me.dv.vcf.gz | grep -E '^SN|^TSTV'

Consider a 30x PCR-free library on GRCh38, where 30x means each position is covered by about thirty independent reads on average. The targets are a mean depth of 28-32x and at least 90% of the callable genome at 20x or better. You also want a genome-wide transition-to-transversion ratio (Ti/Tv) near 2.0-2.1 and 4.4-5.0 million variants. Of those, roughly 3.2-3.6 million should be heterozygous. A Ti/Tv below 1.9 indicates that your filtering let junk through. Coverage is never uniform: GC-rich first exons and high-homology segments drop out systematically. Library preparation, particularly fragmentation chemistry, measurably changes how evenly clinically relevant genes are covered ⁵. Check depth gene by gene before concluding that a gene is clean.

To calibrate the pipeline itself rather than your sample, run the same commands on GIAB HG002 FASTQs. Then compare the result against the benchmark callset using hap.py --engine=vcfeval, restricted to the high-confidence BED. Those reference materials exist precisely so that you can put a number on your own false-positive and false-negative rates instead of guessing ⁶. A healthy 30x pipeline should give SNV F1 above 0.995 and indel F1 around 0.99. If yours is materially worse, fix that before interpreting anything.

Even a well-calibrated short-read genome leaves gaps. Roughly 2-3% of the genome sits in regions where short reads cannot be placed uniquely. Functional interpretation of non-coding variants also remains largely unsolved. The ENCODE pilot showed that most of the genome is transcribed and dense with candidate regulatory elements, which is a statement about how much we can see rather than how much we can interpret ⁷.

9. Handle the file like the identifier it is

Genomic data is unusual among personal records because it implicates people who never handed it over. Your genome identifies your siblings and parents as well as your children and cousins, none of whom consented to anything. Long-range familial search on a database covering about 2% of a target population can return a third-cousin-or-closer match for most individuals of European descent in the United States ⁸. Uploading to a free analysis site is therefore a decision you are making on behalf of your relatives.

Our own practice is straightforward. Keep raw files on an encrypted volume, never upload to a service without a written deletion policy, and run every annotation offline. The broader consent problem is well described in the whole-genome sequencing ethics literature. That literature covers incidental findings and the near-impossibility of meaningful re-consent as variant interpretation changes over years ⁹. Decide in advance which categories of finding you want to see.

Common problems

A handful of failure modes account for most of the trouble people run into, and each has a specific remedy.

Sites silently dropped during --tsv2vcf. This is almost always a build mismatch or a chr prefix mismatch between your TSV and the FASTA. Compare a few known positions against dbSNP by hand to identify which.

A/T and C/G flips. Strand-ambiguous sites cannot be resolved by allele identity alone. When merging with a reference panel, either drop them using plink2 --exclude on a list of ambiguous sites, or resolve them by allele frequency. Never resolve by frequency for variants near 50%, where the two strands are indistinguishable.

Indels encoded as I/D. There is no recovery path for these, since the allele sequence was never recorded. Filter them out.

VEP picking the wrong transcript. The --everything flag emits all consequences. Add --pick_allele_gene --canonical. Better still, use MANE Select transcripts when comparing your results against clinical reports.

ClinVar labels from single submitters with no criteria. A large share of “pathogenic” annotations come from old submissions that carry no assertion criteria. Always carry CLNREVSTAT through your analysis and weight two-star or better.

Missing heritability in polygenic scores. plink2 --score will happily run a European-derived score on any genome. Predictive performance drops substantially outside the ancestry of the discovery GWAS, and the resulting number is not comparable across populations.

Imputation as a substitute for sequencing. Imputing array data against a reference panel gives good accuracy for common variants and poor accuracy for rare ones. Rare variants are the ones you were hoping to see. Imputed rare genotypes should not be interpreted clinically.

Finding something alarming. Stop the analysis, resist the urge to search for treatments, and take the VCF coordinates to a genetic counselor. Confirmation in a clinical lab is the only thing that makes an array or research-pipeline call actionable.

Woolf Software 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. Monkol Lek, Konrad J. Karczewski, Eric Vallabh Minikel, et al. Analysis of protein-coding genetic variation in 60,706 humans. Nature, 2016. https://doi.org/10.1038/nature19057 ↩

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. Swapan Mallick, Heng Li, Mark Lipson, et al. The Simons Genome Diversity Project: 300 genomes from 142 diverse populations. Nature, 2016. https://doi.org/10.1038/nature18964 ↩

4. Maanasa Raghavan, Pontus Skoglund, Kelly E. Graf, et al. Upper Palaeolithic Siberian genome reveals dual ancestry of Native Americans. Nature, 2013. https://doi.org/10.1038/nature12736 ↩

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

6. Justin M. Zook, David Catoe, Jennifer McDaniel, et al. Extensive sequencing of seven human genomes to characterize benchmark reference materials. bioRxiv (Cold Spring Harbor Laboratory), 2015. https://doi.org/10.1101/026468 ↩

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

8. Yaniv Erlich, Tal Shor, Itsik Pe’er, et al. Identity inference of genomic data using long-range familial searches. Science, 2018. https://doi.org/10.1126/science.aau4832 ↩

9. Wim Pinxten, Heidi Howard. Ethical issues raised by whole genome sequencing. Best Practice & Research Clinical Gastroenterology, 2014. https://doi.org/10.1016/j.bpg.2014.02.004 ↩