How to Read Your Raw DNA Data Yourself

Woolf Software
genomics vcf annotation bioinformatics raw-data
A feathered specimen creature on a black plinth, with small glowing markers tagging scattered feathers and large unmarked dark areas.

By the end of this walkthrough you will have converted your consumer genotype file into a coordinate-sorted, reference-checked VCF on GRCh38. It will be annotated with gene names and consequence predictions, plus ClinVar assertions and population allele frequencies. You will also know which of the resulting calls are worth a second look and which are array artifacts.

A VCF, or Variant Call Format file, is the standard text format for recording genetic variants at specific genomic coordinates along with the supporting information about each one. To produce one you need a Linux or macOS machine with about 40 GB of free disk and a GRCh38 reference FASTA. You also need bcftools (1.19+) and samtools, plus either Ensembl VEP or snpEff/SnpSift. Docker is fine and probably easier. You also need your raw data export. 23andMe gives you a genome_<name>_v5_Full_<date>.txt, Ancestry gives you AncestryDNA.txt, and MyHeritage gives you a CSV. All three are tab- or comma-separated tables of roughly 600,000 to 750,000 rows.

One warning belongs before the first command. Nothing below is diagnostic. Consumer arrays are not clinical assays, and they are not validated for the variants they happen to include. A single genotype call is a hypothesis rather than a result. If anything you find looks medically consequential, the path forward is a clinician and a CLIA-certified confirmatory test, not a forum thread.

1. Look at the file before you process it

Before running any tool, spend a minute reading the file itself. The header carries the two facts that determine everything downstream: which genome build the coordinates refer to, and which version of the chip produced them. Open it with head -30.

# This data file generated by 23andMe at: Thu Mar 14 09:12:33 2024
# This file contains raw genotype data, including data that is not used in 23andMe reports.
# Below is a text file containing assembly build 37 (GRCh37) coordinates...
# rsid  chromosome  position  genotype
rs4477212   1   82154   AA
rs3094315   1   752566  AG

Ancestry’s format is similar in content but different in layout. It splits the two alleles into separate columns and uses 0 to mark a position where the array failed to produce a call:

rsid    chromosome  position    allele1 allele2
rs4477212   1   82154   A   A

Three things in this file matter, and it helps to take them one at a time. The first is the build. Almost every consumer export is still on GRCh37, also written hg19, while every modern annotation resource is built on GRCh38. You will therefore have to lift the coordinates over from one build to the other. Liftover is where most home analyses go wrong without announcing it.

The second is strand. 23andMe reports genotypes on the plus strand of the reference assembly, which is convenient. Even so, the mapping from rsID to allele in dbSNP has flipped for a nontrivial number of markers over the years. An rsID alone is therefore not a safe key for matching.

The third is the no-call rate, meaning the number of probes that failed to return a genotype. You can count them directly:

awk '$1 !~ /^#/ && $4 ~ /-|^00$/' genome.txt | wc -l

Typical is 5,000 to 20,000 out of 650,000. If you see 60,000 no-calls, either the sample degraded or the array underperformed. Every downstream conclusion becomes shakier in that case.

2. Convert to VCF against the reference

The next step turns the genotype table into a proper VCF. Resist the temptation to write your own converter. The problem looks like a simple mapping of columns to columns, and it is not. You need to determine which allele is the reference allele at each position and handle hemizygous X and Y calls. You also need to drop the indels that the array encodes as DD/II and flip strands where needed. bcftools convert --tsv2vcf does all of this and checks each call against the actual reference base.

First prepare the input, using the 23andMe layout as the example:

grep -v '^#' genome.txt \
  | awk 'BEGIN{OFS="\t"} $4 !~ /[DI-]/ {print $2,$3,$1,$4}' \
  | sort -k1,1V -k2,2n \
  | bgzip -c > geno.tsv.gz
tabix -s1 -b2 -e2 geno.tsv.gz

Then convert. Note that the reference here is GRCh37, because that is the build the coordinates in the file refer to:

bcftools convert --tsv2vcf geno.tsv.gz \
  -f GRCh37.primary_assembly.fa \
  -s MYSAMPLE \
  -Oz -o raw.b37.vcf.gz
bcftools index raw.b37.vcf.gz

bcftools reports how many sites it skipped because neither reported allele matched the reference base, and that number is a useful sanity check. A few thousand skips out of 650,000 is normal. Those are mostly A/T and C/G markers where the strand is ambiguous, plus positions where the reference itself changed. If the tool skips 100,000 sites, you used the wrong build.

3. Lift over to GRCh38, and check the losses

With a valid GRCh37 VCF in hand, the next task is moving those coordinates onto GRCh38 so that modern annotation resources apply. You can use bcftools +liftover, the plugin from the score collection, or Picard LiftoverVcf with a chain file. Picard is slower, but its rejection log is clearer, and clarity matters more here than speed.

java -jar picard.jar LiftoverVcf \
  I=raw.b37.vcf.gz \
  O=raw.b38.vcf.gz \
  CHAIN=hg19ToHg38.over.chain.gz \
  REJECT=rejected.vcf \
  R=GRCh38.primary_assembly.fa \
  WARN_ON_MISSING_CONTIG=true

Read rejected.vcf rather than skipping past it. You should lose well under 1% of sites. These are concentrated in segmental duplications, the MHC, and a handful of regions that were rearranged between builds. The MismatchedRefAllele rejections are the interesting ones. Those are positions where the reference base itself changed between GRCh37 and GRCh38. Had you lifted over naively by adding a fixed offset, you would now be reading the wrong allele. The MHC on chromosome 6 deserves special suspicion in any personal-genome exercise, because its diversity and the reference bias that follows from it mean array probes and alignments both behave badly there 1.

Finally, sort and index the lifted file, and look at the summary statistics:

bcftools sort -Oz -o final.b38.vcf.gz raw.b38.vcf.gz
bcftools index final.b38.vcf.gz
bcftools stats final.b38.vcf.gz | head -30

4. Annotate

Annotation is the step that turns coordinates and alleles into something interpretable. It tells you which gene a variant falls in and what it does to the protein. It also reports what ClinVar submitters have said about it and how common it is in the population. Ensembl VEP is the tool we would use, because its consequence calls are computed against the same gene models Ensembl publishes. That means you can trace any annotation back to a transcript ID and an exon 2. Run it offline with a cached database, since the online web interface will choke on 600,000 variants.

vep -i final.b38.vcf.gz -o annotated.vcf --vcf \
  --cache --offline --assembly GRCh38 --dir_cache ~/.vep \
  --everything --pick_allele_gene \
  --fasta GRCh38.primary_assembly.fa \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --custom clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
  --fork 4 --compress_output bgzip

Two flags are doing most of the work. --everything pulls in SIFT, PolyPhen, gnomAD frequencies, canonical flags, and HGVS notation. --pick_allele_gene gives one consequence per gene rather than one per transcript, which keeps the output readable. Drop it if you want every transcript.

The annotated file is still far too large to read, so filter it down to something a person can work through. A sensible starting point is ClinVar pathogenic or likely pathogenic calls that carry at least two-star review status. Two stars means multiple submitters agree or an expert panel has reviewed the assertion.

bcftools view -i 'INFO/ClinVar_CLNSIG ~ "athogenic"' annotated.vcf.gz \
  | SnpSift filter "(ClinVar_CLNREVSTAT =~ 'multiple_submitters' | ClinVar_CLNREVSTAT =~ 'reviewed_by_expert_panel')" \
  | bcftools query -f '%CHROM\t%POS\t%ID\t%REF\t%ALT\t[%GT]\t%INFO/ClinVar_CLNSIG\t%INFO/ClinVar_CLNDN\n'

Expect a few dozen hits, and expect nearly all of them to be heterozygous carrier states for recessive conditions. That is exactly what population genetics predicts, since everyone carries a handful of loss-of-function alleles. Treat it as the baseline rather than as a finding.

5. Check genotype, zygosity, and whether the variant is even on the chip

Every hit from the previous step needs to survive four questions before it deserves any attention. Ask them in order, because each one eliminates a different class of false alarm.

The first question is whether the genotype is homozygous reference. A surprising number of “raw data” reports flag a variant when your call is 0/0, because the tool matched on rsID and ignored the GT field entirely. Confirm the actual call directly:

bcftools query -r chr13:32338000-32338200 -f '%POS %REF %ALT [%GT]\n' final.b38.vcf.gz

The second question is whether the condition is recessive and you are heterozygous. Carrier status for a recessive condition is a reproductive-planning fact rather than a personal health prediction. Interpreting it properly is a genetic counselor’s job.

The third question is whether the allele is common. Pull the gnomAD frequency from the VEP output. A variant labeled pathogenic in ClinVar but present at 3% in any population is almost certainly misclassified or of very low penetrance. Frequencies also differ sharply between populations, and most reference panels overrepresent European ancestry. A variant that looks rare globally may therefore be common in your own background 3. Population reference databases such as PGG.Population help you check whether an allele’s frequency is unusual in the ancestry group you belong to rather than in a pooled average 4. Keep in mind too that ancestry is a poor proxy for genotype at any single locus, since variation is overwhelmingly within-group rather than between-group. “This variant is rare in Europeans” tells you very little about you specifically 5.

The fourth question is whether the site is well typed on the array at all. Arrays genotype specific probe positions and do not sequence genes. BRCA1 and BRCA2 together contain thousands of reported pathogenic variants, most of them rare frameshifts and nonsense changes unique to a few families. A consumer chip carries on the order of a few dozen positions in those genes, usually founder variants. A negative result at those positions says almost nothing about the rest of the gene. This is the single largest misreading of raw data we see.

6. Do the things arrays are genuinely good at

Having established what array data cannot support, it is worth being equally clear about what it does well. Arrays are well suited to analyses that aggregate across many markers, because those do not depend on any single call being right.

Ancestry and relatedness are the clearest example. Merge your VCF with 1000 Genomes phase 3, lifted to GRCh38. Prune for linkage disequilibrium so that correlated markers do not dominate, then run a principal components analysis:

plink2 --vcf merged.vcf.gz --set-all-var-ids '@:#:$r:$a' \
  --indep-pairwise 200 50 0.2 --out prune
plink2 --vcf merged.vcf.gz --extract prune.prune.in \
  --pca 10 --out mypca

Plot PC1 against PC2 and you will land where you expect. The reason it works is that the analysis averages over hundreds of thousands of markers, so individual genotyping errors wash out. Deep-ancestry claims about specific haplogroups and migrations belong to a different genre and are much softer than the PCA 6.

Several other analyses share the same statistical advantage. Runs of homozygosity are one. Imputation of untyped common variants is another, using the Michigan or TOPMed imputation server. Imputation returns dosages along with an R² quality measure per variant, and you should discard anything below R² 0.8. Polygenic scores computed from published weights are a third. They are the one output where array data is arguably better than a naive low-coverage genome, because the scores were trained on imputed array data in the first place.

7. Know when to stop and sequence

Eventually the ceiling is the technology rather than the analysis. An array asks about 650,000 pre-chosen yes/no questions. Whole-genome sequencing reads the whole thing. That includes structural variants, everything in genes the chip only samples, and the rare variants that constitute most of the medically interpretable signal. Illumina’s reversible-terminator chemistry established the basic accuracy profile that 30x short-read whole-genome sequencing still runs on. Here 30x means each base is read about thirty times on average, and a modern 30x genome yields roughly 4 to 5 million variants per person instead of 650,000 7. Sequencing carries its own failure modes, including coverage biases tied to GC content and library preparation that vary systematically across the genome 89.

It is also worth being clear about what sequencing does not solve, because interpretation remains the bottleneck. Most variants you find will be of uncertain significance, and the fraction that can be confidently called is small 1. Professional bodies have spent a decade working out how incidental findings from genome sequencing should be handled. The conclusion has consistently been that returning them requires clinical framing and counseling rather than a PDF 10.

Common problems

The failure modes below account for most of the wrong conclusions drawn from consumer genotype data. Each one is easy to check for once you know it exists.

Wrong build, silent failure. This happens when you process a GRCh37 file against a GRCh38 reference. bcftools convert will skip most sites and tell you so, while a hand-rolled script will simply produce nonsense. Always check the site count in bcftools stats against the line count of the input.

Strand flips on A/T and C/G markers. These are the palindromic SNPs, and you cannot resolve them from the genotype alone. If a variant matters, verify it independently. If you are merging with an external dataset, drop them: bcftools view -e 'REF="A" && ALT="T"' etc., or use plink2 —snps-only just-acgt` plus an allele-frequency-based flip check.

rsID drift. dbSNP merges and retires rsIDs over time, so rs1234 in a 2015 export may map to a different merged record today. Match on chromosome and position along with ref and alt, and treat the rsID as a label rather than a key.

No-calls read as reference. A -- or 00 genotype means the array failed at that probe, not that you carry the reference allele. Make sure your conversion emits ./. and not 0/0, and check with bcftools query -f '[%GT]\n' final.b38.vcf.gz | sort | uniq -c.

Trusting an uploaded report. Third-party interpretation sites match rsIDs against literature databases that include underpowered candidate-gene studies from the 2000s, most of which never replicated. A report with 3,000 “findings” has roughly 3,000 false positives in it. Read the raw call and the primary literature, or read nothing.

Mitochondrial and Y confusion. Arrays type a few thousand mtDNA and Y positions, which support haplogroup assignment and little else. Heteroplasmy is invisible to a genotyping array.

Assuming a gene is covered. Before you conclude anything about a gene, count the markers you have in it:

bcftools view -r chr17:43044295-43125483 final.b38.vcf.gz | grep -vc '^#'

Run that on BRCA1 and look at the number. It will calibrate your expectations faster than any argument.

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. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩ ↩2

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

  3. Israel Aguilar-Ordóñez, Eugenio Guzman-Cerezo, David Torres-Treviño, et al. Whole genome sequencing of 1427 Mexican individuals from the oriGen cohort. Nature Communications, 2026. https://doi.org/10.1038/s41467-026-77389-0 ↩

  4. Chao Zhang, Yang Gao, Jiaojiao Liu, et al. PGG.Population: a database for understanding the genomic diversity and genetic ancestry of human populations. Nucleic Acids Research, 2017. https://doi.org/10.1093/nar/gkx1032 ↩

  5. Jeffrey C. Long, Jie Li, Meghan Elisabeth Healy. Human DNA sequences: More variation and less race. American Journal of Physical Anthropology, 2009. https://doi.org/10.1002/ajpa.21011 ↩

  6. Unknown. Deep ancestry: inside the Genographic Project. Choice Reviews Online, 2007. https://doi.org/10.5860/choice.44-5048 ↩

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

  8. Ming-Sin Cheung, Thomas A. Down, Isabel Latorre, et al. Systematic bias in high-throughput sequencing data and its correction by BEADS. Nucleic Acids Research, 2011. https://doi.org/10.1093/nar/gkr425 ↩

  9. Daniel C. Koboldt, Li Ding, Elaine R. Mardis, et al. Challenges of sequencing human genomes. Briefings in Bioinformatics, 2010. https://doi.org/10.1093/bib/bbq016 ↩

  10. Carla van El, Martina C. Cornel, Pascal Borry, et al. Whole-genome sequencing in health care. European Journal of Human Genetics, 2013. https://doi.org/10.1038/ejhg.2013.46 ↩