How to Analyze Your Raw DNA Data Without a Paid Upload Site

Woolf Software
genomics raw-dna vcf annotation privacy
A silver and indigo feathered reptile on a black pedestal, a glowing coiled filament visible beneath its parted breast feathers.

Consumer genotyping files are cheap and portable. They are also far more useful than the trait reports that come bundled with them. This walkthrough shows you how to get that value out yourself.

By the end you will have converted your consumer array file into a sorted, reference-checked VCF on GRCh38. A VCF, or Variant Call Format file, is the standard way genomics tools represent a person’s variants against a reference genome. You will have annotated that file with ClinVar and Ensembl VEP, scored it against a few PGS Catalog models, and optionally imputed it to roughly 40M sites on a free server. You will also have a clear picture of which upload sites are worth your file and which are a data-collection funnel with a trait report bolted on.

The prerequisites are modest. You need three things:

  • Your raw download from 23andMe or AncestryDNA. MyHeritage, FamilyTreeDNA, and Living DNA files also work. It arrives as a 15-25 MB text file of roughly 600k-700k lines.
  • A Linux or macOS machine with about 50 GB free.
  • bcftools and plink2 installed. You also need plink 1.9, CrossMap, and Docker.

Everything below is measurement and interpretation. Nothing here is medical advice. Any variant you would act on needs confirmation in a clinical lab and a conversation with a genetics professional.

1. Read your file before you upload it anywhere

The first task is simply to understand what you have. Ten minutes of inspection before you hand the file to any website will save you from a great deal of confusion later.

Every vendor ships a slightly different tab-delimited or CSV format. As of this writing they are all on GRCh37, also called build 37 or hg19. That is the human reference assembly released in 2009. They all report genotypes on the plus strand of that reference.

A 23andMe file from the v5 chip covers roughly 630k sites and looks like this:

# rsid	chromosome	position	genotype
rs548049170	1	69869	TT
rs13328684	1	74792	--

AncestryDNA covers about 670k sites and splits the genotype across two allele columns. It also uses numeric chromosome codes, in which 23=X, 24=Y, 25=PAR, and 26=MT:

rsid	chromosome	position	allele1	allele2
rs4477212	1	82154	A	A

The others are variations on the same theme. MyHeritage is a quoted CSV with the columns RSID,CHROMOSOME,POSITION,RESULT, and FamilyTreeDNA is similar. Living DNA is tab-delimited with a single genotype column.

Whatever the vendor, the same first pass applies. It is worth running before any website sees the data:

grep -vc '^#' genome.txt                      # ~630000 for 23andMe v5
awk '!/^#/ && $4=="--"' genome.txt | wc -l    # no-calls; expect 0.3-2%
awk '!/^#/{print $2}' genome.txt | sort | uniq -c   # chromosome coverage
sha256sum genome.txt > genome.sha256

Those four commands tell you the site count and the fraction of positions the chip failed to call. They also give you the per-chromosome coverage and a checksum. The checksum lets you prove later that you are working with the original download. A no-call rate above about 3% means the array underperformed, and every downstream analysis inherits the gaps.

Two more cheap sanity checks are worth doing. The first is X heterozygosity. A male sample should show almost no heterozygous calls of the AG variety on chromosome X outside the pseudoautosomal regions, the small stretches that chromosome X shares with the Y. The second is overall heterozygosity, which sits near 30-35% of called autosomal sites for most samples. A wildly off value usually means a mangled download, or a file that has somehow been merged from two people.

2. Pick upload sites deliberately, and know what each does

There is no shortage of sites that will take your raw file for free. They are doing quite different things with it, and the trade you are making differs accordingly. The free sites cluster into three categories.

The first category is genealogy matching. GEDmatch is the one that matters here. You upload a raw file, it produces a kit number, and you can then run one-to-many matching against everyone else who has uploaded. It also offers segment-level one-to-one comparisons and admixture models. The core tools are genuinely free. Tier 1, which adds triangulation and some visualization, is a paid subscription. GEDmatch is owned by Verogen/QIAGEN, and every kit carries a law-enforcement matching flag that you set at upload.

Other matching services take uploads on similar terms. FamilyTreeDNA accepts free uploads and gives you the match list, charging to unlock ethnicity estimates and the chromosome browser. MyHeritage accepts uploads and provides matching plus ethnicity at no cost. Living DNA accepts uploads and returns a free ancestry report. AncestryDNA and 23andMe accept no uploads at all, in either direction. So the answer to “can I upload my raw DNA to 23andMe” is no. You have to test with them.

The second category is interpretation. Genetic Genie runs free methylation and detox panels off a handful of rsIDs, the familiar MTHFR and COMT set among them. It now offers a broader “Discovery” variant report as well. Promethease was the serious entry in this category, cross-referencing your genotypes against SNPedia. It costs a few dollars, and its future has been uncertain since the MyHeritage acquisition. Codegen, Genomelink, and Xcode Life hand you a free teaser report and then upsell additional panels. Similar sites work the same way. None of these tools is a substitute for running annotation yourself, and step 5 reproduces most of what they do with much better provenance.

The third category is research and open data. openSNP lets you publish your file into the public domain alongside phenotype answers. That is a different bargain from the others, because you are contributing to a participant-driven research commons rather than buying a report 1. The ethics literature on these projects is worth reading before you post an irrevocably public genome 1.

One consideration cuts across all three categories. Your file is an identifier for you and a partial identifier for every first- and second-degree relative you have, none of whom consented to the upload. Terms of service on consumer genomics sites change with ownership. The legal protections for direct-to-consumer genetic data are thinner and more fragmented than most people assume 2. Users also tend to treat these files as living, re-interpretable objects that they return to for years, which means the upload decision is not a one-time disclosure 3. Our own practice is to upload to GEDmatch when genealogy is the goal, and otherwise to keep the file local.

3. Convert to VCF against a real reference

Text genotypes are perfectly adequate for match sites and useless for annotation. The next step is therefore a single careful conversion to VCF, the standard Variant Call Format used by essentially every genomics tool.

The important part of the conversion is checking each call against the reference base at that position. Vendor files contain sites where the listed allele does not match GRCh37. The causes include strand ambiguity, indel representation, and positions that have moved between assembly versions.

Start by obtaining a reference FASTA that matches your file’s build. Either human_g1k_v37.fasta or the GRCh37 primary assembly will do. Then run:

bgzip -c genome.txt > genome.txt.gz
bcftools convert --tsv2vcf genome.txt.gz \
  -f human_g1k_v37.fasta -s ME \
  -c ID,CHROM,POS,AA \
  -Oz -o me.b37.vcf.gz
bcftools index -t me.b37.vcf.gz

bcftools convert --tsv2vcf expects the 23andMe layout by default. If you are working from an AncestryDNA file, collapse the two allele columns and translate the numeric chromosome codes first:

awk 'BEGIN{OFS="\t"} !/^#/ && $1!="rsid" {
  c=$2; if(c==23)c="X"; else if(c==24)c="Y"; else if(c==25)c="X"; else if(c==26)c="MT";
  print $1, c, $3, $4 $5 }' AncestryDNA.txt | bgzip > ancestry.tsv.gz

Read the conversion log, which reports how many sites were skipped for reference mismatch. A few thousand out of 650k is normal. Tens of thousands means you used the wrong build or the wrong chromosome naming, and you should fix that before going further.

If you would rather stay in PLINK-land for scoring and kinship work, the equivalent path is:

plink --23file genome.txt ME ME i --output-chr MT \
      --snps-only just-acgt --make-bed --out me_b37
plink2 --bfile me_b37 --fa human_g1k_v37.fasta --ref-from-fa force \
       --make-pgen --out me_b37_ref

The --ref-from-fa force flag is what stops PLINK from guessing the reference allele by allele frequency. Left to guess from a single sample, it silently flips alleles and quietly corrupts every polygenic score you compute afterward.

4. Lift to GRCh38

Almost every current annotation resource is built GRCh38-first, so your coordinates need to move to that assembly. CrossMap does this using the UCSC chain file, which maps positions between the two builds:

CrossMap vcf hg19ToHg38.over.chain.gz me.b37.vcf.gz \
  GRCh38.primary_assembly.genome.fa me.b38.vcf
bcftools sort me.b38.vcf -Oz -o me.b38.vcf.gz && bcftools index -t me.b38.vcf.gz
bcftools norm -c x -f GRCh38.primary_assembly.genome.fa -m -any \
  me.b38.vcf.gz -Oz -o me.b38.norm.vcf.gz

Expect to lose roughly 0.2-0.5% of sites to regions that do not map between builds. Expect CrossMap to emit unsorted output, which is why the explicit bcftools sort is there. The -c x argument to norm rechecks the reference allele against the new build and drops sites that disagree.

Keep the unmapped file rather than discarding it. If a site you care about turns up in it, look the site up by rsID in dbSNP before assuming it is gone.

5. Annotate with VEP and ClinVar

This is the step that replaces Promethease. It carries the considerable advantage that you know exactly which database version produced each line of output.

Run Ensembl’s Variant Effect Predictor in a container. That way the tool, its cache, and its plugins are all pinned to known versions 4:

docker run -v $PWD:/data ensemblorg/ensembl-vep \
  vep --offline --cache --dir_cache /data/vep_cache \
      --assembly GRCh38 --everything --vcf --compress_output bgzip \
      --custom /data/clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
      -i /data/me.b38.norm.vcf.gz -o /data/me.vep.vcf.gz

Pull clinvar.vcf.gz from the NCBI FTP site and record its release date alongside your results. ClinVar classifications change over time, so the version matters. Then extract anything ClinVar calls pathogenic with a meaningful review status:

bcftools +split-vep me.vep.vcf.gz -f '%CHROM %POS %ID %REF %ALT [%GT] %SYMBOL %ClinVar_CLNSIG %ClinVar_CLNREVSTAT\n' \
  -i 'ClinVar_CLNSIG ~ "athogenic" && ClinVar_CLNREVSTAT ~ "multiple_submitters"' \
  -d | column -t

On a 650k-site array you will typically get somewhere between a handful and a few dozen hits. Most will be heterozygous carrier states for recessive conditions.

Two caveats matter more than the results themselves. The first concerns the technology. Consumer arrays are optimized for common variation, and the rare pathogenic variants they do probe are exactly where genotyping calls are least reliable. Cluster separation for a variant present in 1 in 20,000 people is poor. A pathogenic call from an array file is a hypothesis, and it should be treated as one.

The second caveat concerns the process. Clinical interpretation is a curated, versioned, laboratory-controlled process with confirmation and reporting standards that a local VEP run does not replicate 5. Access to your own sequence data is genuinely valuable. The research literature supports giving participants that access while also documenting the misinterpretation risk that comes with it 6.

If an annotation suggests a pathogenic variant with any clinical weight, treat it as unconfirmed until a CLIA or equivalent laboratory retests the site and a clinical geneticist or genetic counselor interprets it.

Pharmacogenomics is a common reason people upload files in the first place, and it can be handled locally too. PharmCAT will take your VCF and produce star-allele calls for CYP2C19, CYP2D6, SLCO1B1, and others. The CYP2D6 result is only partial, since arrays miss the structural variation in that gene. Read the output as a description of your genotype and nothing more. Any prescribing implication is a conversation with a physician or clinical pharmacist, full stop.

6. Polygenic scores worth computing

Polygenic scores are easy to compute and easy to compute wrongly, so the setup matters. Download harmonized scoring files from the PGS Catalog. They are named in the pattern PGS000xxx_hmPOS_GRCh38.txt.gz and give you rsID, effect allele, and weight already mapped to build 38.

plink2 --pfile me_b38 \
  --score PGS000018_hmPOS_GRCh38.txt 1 4 6 header-read \
          cols=+scoresums,+denom no-mean-imputation \
  --out pgs000018

Column indices differ from file to file, so check the header before running. Use no-mean-imputation and then look at the denominator that --score reports. Suppose the model uses 1.1M variants and your array covers 180k of them. Your score is then not comparable to the published distribution. In that case you should either impute first (step 7) or choose a score built on a smaller variant set.

Scores are also calibrated within a specific ancestry group and shift badly outside it. An absolute percentile means little unless your genetic ancestry matches the training cohort. We treat single-person PGS output as a relative ranking within a model, nothing more.

7. Impute, if you want the extra sites

Imputation fills in genotypes at sites your chip never measured, using patterns of linkage from a large reference panel. The TOPMed Imputation Server is free with registration and takes per-chromosome bgzipped GRCh38 VCFs. It returns roughly 40-50M imputed sites with an R2 quality field per variant.

Split your file and check it before submitting:

for c in {1..22}; do
  bcftools view -r chr${c} me.b38.norm.vcf.gz -Oz -o chr${c}.vcf.gz
done

Run the server’s pre-check script, either the HRC/1000G checking tool or checkVCF.py. It catches allele-frequency outliers and strand problems before submission rather than after a failed run.

Afterward, filter the output hard. Use bcftools view -i 'R2>0.8' for common variants, and expect imputation of low-frequency variants to be poor regardless of the threshold you pick. Imputation gives you genotype probabilities at sites nobody measured on you. That is useful for polygenic scoring and worthless for anything that hinges on a single rare variant.

8. Do your own relative matching

If you have raw files from several family members, you do not need a match site to compute how closely they are related. Merge the files down to their intersecting sites and run:

plink2 --pfile merged --make-king-table --out kin

KING kinship coefficients land near 0.25 for parent-child and full-sibling pairs. Those two are distinguished by the IBS0 column, which sits near zero for parent-child. Half-siblings, grandparent-grandchild, and aunt or uncle relationships come in near 0.125. First cousins land near 0.0625. This is exactly the calculation behind GEDmatch’s one-to-one tool, minus the shared database. That database is the only reason to upload at all.

Common problems

Most failures in this pipeline are coordinate or allele bookkeeping errors, and they tend to be silent rather than loud. The following are the ones we see most often, along with how to recognize and fix each.

  • Reference mismatch on thousands of sites. You converted a build 37 file against a build 38 FASTA, or the reverse. Check a known position: rs1815739 sits at chr11:66,560,624 in GRCh37 and chr11:66,793,158 in GRCh38.
  • Chromosome naming mismatch. A file using 1 where a tool expects chr1 breaks bcftools annotate, CrossMap, and the imputation server in different ways. Fix it once with bcftools annotate --rename-chrs.
  • Allele flipping after PLINK conversion. Without --ref-from-fa force, PLINK sets the reference allele from your single sample. Every homozygous-alternate site then becomes the “reference” and your scores invert. Verify by spot-checking three known sites against dbSNP.
  • Indels and CNVs missing entirely. Arrays genotype a fixed set of mostly single-nucleotide sites. Arrays genotype a fixed set of mostly single-nucleotide sites. Structural variation, repeat expansions, and most indels are absent from the file. Their absence is not a negative result. Short-read whole-genome sequencing recovers classes of variation that targeted assays never see 7, which is a large part of why a genotyping file cannot be treated as a genome.
  • No-calls read as homozygous reference. A -- in the source file should become a missing genotype rather than 0/0. Confirm that bcftools stats reports a nonzero missing count.
  • A surprising or alarming result. Array-based calls for rare variants carry a meaningful false-positive rate. Quality assurance procedures exist in clinical sequencing precisely because unvalidated pipelines produce confident wrong answers 8. Confirm in a clinical lab before you believe it, and take the result to a clinician rather than to a forum.
  • Uploading to sites you cannot later exit. Deletion policies vary and ownership changes. Decide once, per site, with the understanding that a match database is only useful because other people made the same disclosure 2.

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. Michelle L. McGowan, Suparna Choudhury, Eric T. Juengst, et al. “Let’s pull these technologies out of the ivory tower”: The politics, ethos, and ironies of participant-driven genomic research. BioSocieties, 2017. https://doi.org/10.1057/s41292-017-0043-6 ↩ ↩2

  2. Henry T. Greely. The Future of DTC Genomics and the Law. The Journal of Law Medicine & Ethics, 2020. https://doi.org/10.1177/1073110520917003 ↩ ↩2

  3. Minna Ruckenstein. Keeping data alive: talking DTC genetic testing. Information Communication & Society, 2016. https://doi.org/10.1080/1369118x.2016.1203975 ↩

  4. Konstantinos Krampis, Tim Booth, Brad Chapman, et al. Cloud BioLinux: pre-configured and on-demand bioinformatics computing for the genomics community. BMC Bioinformatics, 2012. https://doi.org/10.1186/1471-2105-13-42 ↩

  5. Hana Zouk, Eric Venner, Niall J. Lennon, et al. Harmonizing Clinical Sequencing and Interpretation for the eMERGE III Network. The American Journal of Human Genetics, 2019. https://doi.org/10.1016/j.ajhg.2019.07.018 ↩

  6. Susanne B. Haga, Bethany Friedman, Gabriele Richard. Considering the Benefits and Risks of Research Participants’ Access to Sequence Data. Genetic Testing and Molecular Biomarkers, 2017. https://doi.org/10.1089/gtmb.2017.0143 ↩

  7. Australian Pancreatic Cancer Genome Initiative, Nicola Waddell, Marina Pajic, et al. Whole genomes redefine the mutational landscape of pancreatic cancer. Nature, 2015. https://doi.org/10.1038/nature14169 ↩

  8. Laila Sara Arroyo Mühr, Daniel Guerendiain, Kate Cuschieri, et al. Human Papillomavirus Detection by Whole-Genome Next-Generation Sequencing: Importance of Validation and Quality Assurance Procedures. Viruses, 2021. https://doi.org/10.3390/v13071323 ↩