Skip to content

Genetic software: what to run on your own genome

Oak
An iridescent feathered reptilian creature on a black plinth in a dark studio, with bands of blue light sorting rows of its feathers.

Yes, you can sequence and analyze your own genome without a physician in the loop, and the software to do it is open source, well benchmarked, and runnable on a single machine with 64 GB of RAM and a few terabytes of disk. The stack we would run on a 30x whole-genome sample is below.

  • bwa-mem2 for alignment
  • DeepVariant for small variants
  • Manta or dysgu for structural variants
  • ExpansionHunter for repeats
  • VEP for annotation
  • PharmCAT for star alleles
  • pgsc_calc for polygenic scores

Turning a variant call into a clinical decision requires a clinician and, for anything with reproductive or diagnostic weight, a genetic counselor. You cannot do that step alone.

Start by knowing which file you have

Most people asking about genetic software have one of three artifacts, and they are not interchangeable.

A consumer array export is a tab-delimited text file with roughly 600,000 to 900,000 genotyped sites, each a rsID plus an A/B call. It is about 25 MB. It captures common variation well, because the chips are designed around common SNP maps of the kind first assembled two decades ago 1. It is bad at rare variation, which is the part that usually matters for a single individual, and it does not observe indels, structural variants, or anything not on the chip.

A 30x whole-genome sample gives you FASTQ (roughly 100 GB gzipped, paired-end), a CRAM (15–20 GB against GRCh38), and a VCF or gVCF with on the order of 4–5 million small variant calls per individual. This is what you want if you plan to do your own analysis, because every downstream tool can be re-run when the reference, the caller, or the variant database changes.

An exome sits between the two: about 35–50 Mb of target, 6 GB CRAM, ~25,000 coding variants. Good for coding interpretation, blind to everything regulatory, and regulatory variation is most of the genome’s functional content 2.

Ask your provider for FASTQ or CRAM, in writing, before you pay. Some labs deliver only a filtered VCF, and a filtered VCF is a dead end: you cannot recall, you cannot recover reference-confident blocks, and you cannot tell a missing call from a homozygous-reference call.

The stack

Align to GRCh38 with alt contigs and the HLA decoys. Do not align to hg19. The alt-aware path matters for the MHC and for a handful of segmental duplications.

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:s1\tSM:you\tPL:ILLUMINA\tLB:lib1' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa R1.fq.gz R2.fq.gz \
| samtools sort -@ 8 -m 2G -o you.bam -
samtools index -@ 8 you.bam
mosdepth -n --fast-mode --by 1000 you you.bam

Check mosdepth output before doing anything else. You want mean autosomal depth above 28x and at least 90% of callable bases at ≥20x. Below 20x, heterozygous calls start dropping out, and a missed het reads as a reassuring homozygous reference.

For small variants, DeepVariant beats GATK HaplotypeCaller on single-sample Illumina WGS, mostly on indels in homopolymers:

run_deepvariant \
  --model_type=WGS \
  --ref=GRCh38.fa \
  --reads=you.bam \
  --output_vcf=you.vcf.gz \
  --output_gvcf=you.g.vcf.gz \
  --num_shards=32

Then verify. Do not skip this. Run a Genome in a Bottle sample through the identical pipeline (HG002 is the usual choice) and score with hap.py against the v4.2.1 truth set restricted to the high-confidence BED. If your SNV F1 is not above 0.995 and your indel F1 above 0.99, something in your pipeline is wrong, and you will not find out any other way.

Structural variants need a separate caller. Manta (configManta.py --bam you.bam --referenceFasta GRCh38.fa) is fast and conservative on deletions and duplications. dysgu picks up more insertions. Short reads miss large insertions and anything inside repeats, so treat a negative SV result as uninformative rather than reassuring. Repeat expansions need ExpansionHunter with the published variant catalog, and even then the long alleles at loci like FMR1 and C9orf72 are underestimated by short reads.

Annotation is where most self-analysis goes wrong. Use VEP offline with a local cache:

vep -i you.vcf.gz -o you.vep.vcf --vcf --cache --offline \
  --assembly GRCh38 --everything --fork 8 \
  --custom clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.snv.hg38.vcf.gz

Filter on CLNREVSTAT as well as CLNSIG. Filtering on CLNSIG alone is not enough. A “pathogenic” call with one submitter and no assertion criteria is a different object from a reviewed expert panel classification. ClinVar reclassifies variants continuously, and a variant of uncertain significance today may be benign in eighteen months. Re-run annotation quarterly against a fresh ClinVar release. That is the single highest-value recurring job in a personal pipeline.

Pharmacogenomics is the one domain where the software gives you something usable and well-structured. PharmCAT takes a preprocessed VCF and emits CPIC-aligned diplotypes for a set of genes. The set covers CYP2C19, CYP2C9, DPYD, TPMT and SLCO1B1, among others. The known hole is CYP2D6, which has structural rearrangements and hybrid alleles that short reads call poorly. Use PyPGx or the StellarPGx pipeline for CYP2D6 specifically and expect ambiguity. Do not change any medication based on this output. Take the report to a prescriber.

Polygenic scores: run pgsc_calc (Nextflow) against PGS Catalog scoring files rather than hand-rolling a sum. It handles strand flips, allele matching, and ancestry projection against 1000 Genomes reference PCs. The results are population-relative percentiles whose calibration degrades outside the ancestry of the discovery GWAS. Large multi-ancestry sequencing efforts have made this explicit: variant effects on quantitative blood traits differ in frequency and effect across populations, and scores built in one group transfer imperfectly to another 3.

The failure modes that will bite you

Chromosome naming. chr1 versus 1 will silently produce zero annotations. Check with bcftools view -h before blaming the tool.

The MHC, the KIR locus, and the immunoglobulin loci are not reliably callable from short reads aligned to a single linear reference. If you care about HLA typing, run a dedicated typer (HLA-LA or OptiType) rather than reading genotypes off the VCF.

Mitochondrial variants need their own workflow. Heteroplasmy fraction is the quantity of interest, and a diploid caller will not report it correctly. Use the Mutect2 mitochondria mode.

Non-coding variants. You will have millions of them, essentially all uninterpretable at individual level. The ENCODE pilot showed how much of the genome is transcribed or bound by regulatory proteins without a clear mapping to phenotype 2. CADD and SpliceAI give you a prioritization score. They do not give you a verdict.

Incidental findings. A whole genome will show you carrier status for recessive conditions, occasional high-penetrance cancer or cardiac variants, and results with implications for relatives who did not consent to anything. Such findings are a well-described consequence of sequencing rather than targeted testing, and the ethics literature on them predates consumer availability 4. Decide in advance which categories you want to look at, and write the decision down before you run the annotation. The framework for pre-test decision-making and counseling was laid out in detail decades ago and still holds 5.

Where a clinician is not optional

Three cases. First, any variant you intend to act on medically: a confirmatory clinically validated assay in an accredited lab is required, because research-grade calls have a nonzero false positive rate and no chain of custody. Second, reproductive decisions. Carrier screening by NGS across couples is a mature clinical workflow with defined gene panels and interpretation standards. The output is a joint risk assessment. It is not a single genotype 6. Third, suspected Mendelian disease in yourself or a family member. Diagnostic yield depends heavily on phenotype-driven gene selection, and the interpretation of a specific variant in a specific gene often hinges on clinical examination findings that no pipeline sees 7.

Questions people also ask

Can I do genetic testing on my own? You can order sequencing and analyze the output yourself with open-source tools. The analysis is genuinely doable by a competent engineer. The interpretation of a specific finding for your health is not, and confirmatory clinical testing is a separate step.

Can I get a genetic test without a doctor? In most US states, direct-to-consumer sequencing is available without a physician order, sometimes with a review physician the lab supplies. The resulting data is research-grade unless the lab states otherwise. Clinical action on it requires an accredited confirmatory test.

What is the best at-home genetic test? For anyone who wants to do their own analysis: 30x whole-genome sequencing with FASTQ or CRAM delivery. Arrays are cheaper and cover common variation reasonably, but they interrogate under a million fixed sites chosen from common SNP maps 1 and cannot see the rare variation that dominates individual-level findings.

Do I need GRCh38 or T2T-CHM13? Use GRCh38 as the primary, because every annotation database is keyed to it. Run T2T-CHM13 as a secondary alignment if you care about centromeric, acrocentric, or segmentally duplicated regions, which reference-quality work on individual chromosomes showed are the hardest parts to resolve 8.

How much compute does this take? Alignment of a 30x genome is roughly 4–6 hours on 32 cores. DeepVariant is 2–4 hours on the same hardware, under an hour on a single modern GPU. Budget 2 TB of scratch disk and keep the CRAM. Do not keep the BAM.

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

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

  3. Anna V. Mikhaylova, Caitlin McHugh, Linda M. Polfus, et al. Whole-genome sequencing in diverse subjects identifies genetic correlates of leukocyte traits: The NHLBI TOPMed program. The American Journal of Human Genetics, 2021. https://doi.org/10.1016/j.ajhg.2021.08.007 ↩

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

  5. Institute of Medicine. Assessing Genetic Risks. National Academies Press eBooks, 1994. https://doi.org/10.17226/2057 ↩

  6. Julio Martín, Asan, Yuting Yi, et al. Comprehensive carrier genetic test using next-generation deoxyribonucleic acid sequencing in infertile couples wishing to conceive through assisted reproductive technology. Fertility and Sterility, 2015. https://doi.org/10.1016/j.fertnstert.2015.07.1166 ↩

  7. Edwin M. Stone. Genetic Testing for Inherited Eye Disease. Archives of Ophthalmology, 2007. https://doi.org/10.1001/archopht.125.2.205 ↩

  8. Masahira Hattori, Asao Fujiyama, Todd D. Taylor, et al. The DNA sequence of human chromosome 21. Nature, 2000. https://doi.org/10.1038/35012518 ↩

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, blood RNA sequencing, and plasma proteomics, delivered as a detailed report plus every data file. Build your profile.