Reading My Own Genome

A few years ago I sent a tube of saliva to Nebula Genomics. They sent back four files. A 50-something-gigabyte CRAM (the alignment of every read to the reference genome) plus its .crai index, and a 255 megabyte VCF (the variants they called from those reads) plus its .tbi index. The website wraps it in a friendly health summary covering, what, fifteen traits? Twenty? But the actual data, the part I paid for, is those four files. They are objectively interesting and almost completely unreadable to a non-bioinformatician.

The CRAM is the long version of the answer. Every read off the sequencer, where it landed in the genome, what the base call was at every position. You need it if you want to recheck a specific variant, or run a different caller, or look at coverage. For ordinary "what's in my genome" questions it's overkill, and reading it requires samtools and a reference file the same size.

The VCF is the short version. About 5 million variants. Each line says: at this position on this chromosome, the reference is this base, your genotype is this. That's it. No interpretation, no clinical context, no "by the way, this one matters." That part is supposed to come from cross-referencing the VCF against an evidence database, and that's the part Nebula doesn't really do.

So I built it. The tool runs locally, never uploads anything, and pulls in ClinVar for clinical significance, PharmGKB for drug response, gnomAD for population frequencies, PharmCAT for pharmacogene diplotype calling, MyVariant.info for live functional scores, the PubMed literature index, and ten polygenic risk scores from the PGS Catalog. It also does an offspring inheritance simulator, validated against the Genome in a Bottle Consortium's published trio data. This article is what I learned along the way and a few things about my own genome that I did not know last week.

The first finding that mattered

The first thing the tool surfaced that actually changed something for me was a single SNP. Chromosome 10, position 94781859. rs4244285. Heterozygous. The gene is CYP2C19, a cytochrome P450 enzyme that handles roughly 10 percent of all prescription drugs.

CYP2C19 has a star-allele system. *1 is the reference, *2 is a loss-of-function allele caused by exactly the variant I have, *17 is a gain-of-function allele in the promoter region. One *2 means I'm *1/*2, an Intermediate Metabolizer. PharmCAT (the gold-standard pharmacogene caller from Stanford and Penn) confirmed it independently when I ran it locally. CPIC publishes guidelines for what an Intermediate Metabolizer should do for a long list of drugs.

If I'm ever prescribed clopidogrel after a cardiac event, the prodrug doesn't get converted to its active form efficiently, and standard antiplatelet therapy is documented to be less protective in carriers like me. If I'm prescribed citalopram or escitalopram for depression, the drug accumulates higher than expected and the standard starting dose is probably too high for me. Same goes for amitriptyline, sertraline, and proton pump inhibitors at chronic doses.

This isn't speculative. It's a CPIC Level A guideline for my exact diplotype. I've been to a half dozen doctors in the last decade and none of them ever ordered a pharmacogenomic panel, which is the thing that would have surfaced this in any clinical setting. The information was in the VCF the whole time. The bottleneck was the tool to read it.

The Punnett square is correct, mostly

The next layer down is inheritance. Most people remember Punnett squares from ninth grade, draw a 2x2 grid, two parents two alleles each, four boxes for the child. The intuition is right. The math is one step simpler than the grid suggests.

Each parent transmits one of their two alleles to the child with 50/50 probability. The chance the child gets allele A from a given parent is just that parent's dosage divided by two: zero, one half, or one. The joint probability of any specific child genotype is the product of two independent Bernoulli draws. That's it.

You probably learned this with the four-square grid because the grid is easier to draw than the formula. Both give the same answer. Two heterozygous parents produce an affected child 25 percent of the time.

What the textbook version skips is everything that happens once you have multiple genes at once, with linkage between nearby variants, with genotype-quality uncertainty in a real VCF, and with the X-chromosome wrinkle for X-linked conditions. The simulator I built does the textbook math exactly for single variants and switches to Monte Carlo only where it actually adds value: polygenic risk score distributions, multi-locus trait predictions, joint probabilities over many sites at once.

When the answer changes between snapshots

Here's the thing nobody warned me about. ClinVar updates monthly. Variants get reclassified as new evidence comes in, as different labs submit interpretations, as expert panels review old calls. It's a living database, and the pace of change is faster than most people building consumer tools account for.

The first time I ran my tool, the headline finding was a homozygous Pathogenic variant in LDLR, the gene for the LDL receptor, classification "Familial hypercholesterolemia." This is a real thing with a real management pathway. I sat with that for about ten minutes before doing what I should have done first, which was check the live database. The variant had been reclassified to Benign by multiple submitters in the months since I'd last refreshed my local copy.

It wasn't just LDLR. DHCR7 had flipped (Smith-Lemli-Opitz). KCNH2 had flipped (Long QT). FAM161A had flipped (retinitis pigmentosa). Four of my "top findings" were no longer findings at all. They were just old database entries.

Now the tool checks ClinVar live for every actionable variant on every scan. The cached database is for performance. The live query is for trust. Anyone who ships a consumer-facing genomics product with a static interpretation database is shipping answers that are stale within a few months, and I don't think most users realize this.

Validating the simulator against a real family

The offspring simulator is the part I worried about most. The math is straightforward but math compiling clean is not the same thing as math predicting reality. Two hundred and eighteen unit tests across thirteen suites is necessary. It is not sufficient.

The Genome in a Bottle Consortium publishes whole-genome benchmark VCFs for several family trios. The Ashkenazi trio is HG002 (the son), HG003 (the father), HG004 (the mother), and it's the most-cited reference set in the field. I downloaded all three. Then I ran the simulator on the parents only, and compared its predicted distribution at each of HG002's four million non-reference variants against his actual observed genotype.

The metric I cared about most was the average predicted probability the simulator assigned to what was actually observed. If you predicted uniformly across , you'd average 33 percent. A perfectly calibrated simulator would assign high probability to actual observations on average. The Ashkenazi trio came back at 61.0 percent. The independent Han Chinese trio (HG005, HG006, HG007) came back at 62.6 percent. Roughly twice the uniform baseline.

The conditional Mendelian sanity check came out the way Mendel said it should. When both parents are heterozygous at a position and we condition on the child being non-reference (which is what the VCFs actually contain), 33.3 percent of those should be homozygous-alt and 66.7 percent should be heterozygous. Ashkenazi observed: 32.7 percent. Chinese observed: 32.3 percent. Within a percentage point. The small gap is consistent with the known heterozygote calling bias in consensus VCFs, where borderline het calls occasionally round up.

This is what validation looks like when there's a real ground truth available. Here we got a multi-generation family with publicly released data and the simulator's predictions were demonstrably consistent with how that family's genome actually inherited.

How the Monte Carlo holds up

Polygenic risk scores can't be done with exact math. The score is a sum across thousands of variants, each contributing a small weighted dose, and the joint distribution lives in a space too high-dimensional to enumerate. Monte Carlo is the right tool. The question is whether it converges at the rate it's supposed to.

For an unbiased sampler, the deviation from the analytical mean should fall as one over the square root of N. Doubling N cuts error by roughly 1.41. Multiplying by ten cuts by roughly 3.16. I tested the simulator against 10,000 parental loci from the trio data. At N=1,000, mean absolute deviation between empirical and analytical genotype probabilities was 0.0099. At N=10,000 it was 0.0031. At N=100,000 it was 0.0010. The reduction factors were 3.16 and 3.09, basically textbook. KL divergence dropped from 4.3e-4 to 4.3e-6 across the same range. Wilson score 95 percent intervals contained the truth 94.7 to 95.2 percent of the time across 100 independent trials over 500 non-boundary loci.

This is the part that's boring to talk about and important to get right. Every polygenic risk number the tool reports has an error bar attached, and the error bar is one I can vouch for.

What I actually learned about myself

The tool produces 1,073 ClinVar matches across 773 genes from my VCF. After live cross-checking, the actionable set comes down to seven PharmCAT diplotype calls, twelve PharmGKB Level 1A or 1B annotations matching my exact genotype, and a handful of lifestyle-relevant variants the rest of the medical system has never bothered to mention.

I'm a homozygous T at rs5751876 in ADORA2A, which is the variant most associated with caffeine-induced anxiety. I had always blamed mild anxiety on the rest of my brain. Apparently it's also partly that. I'm a heterozygous carrier of ADH1B*2 (rs1229984 CT), the East-Asian-typical fast-metabolism allele, which makes alcohol metabolism roughly 40 to 100 times faster than reference and probably accounts for a flushing reaction I have always sort of ignored. I have one MTHFR C677T allele (mild folate-metabolism reduction; the consensus clinical recommendation is to do nothing). My CYP1A2 puts me at intermediate caffeine clearance. None of this is dramatic individually. Cumulatively it's a coherent picture of why I metabolize substances the way I do.

Eye color predicted brown, 85 percent confidence, off HERC2 alone. Hair predicted dark brown to black. Cilantro tastes fine to me (rs72921001 CC, no soap perception). I can smell asparagus in pee (rs4481887 heterozygous, one functional copy of OR2M7).

What this changes

A whole-genome VCF is decision-relevant data sitting on your hard drive that you can't read. The bottleneck has never been the data or the science. Open-source clinical-grade variant callers exist. Evidence databases are public. Polygenic risk scores are downloadable. The bottleneck is the integration. There's no consumer product that takes your raw VCF and gives you a clinical-grade interpretation, because the regulatory shape of that product would be FDA Software-as-a-Medical-Device, and nobody wants to ship that.

The version I built is for me. It runs on my laptop. My genome never leaves the machine. When a variant gets reclassified, I find out the next time I run a scan. When a doctor mentions a drug I've never taken, I can look up my pharmacogenomics in two minutes. When ClinVar updates, my answers update.

The thing I keep coming back to as the next direction is integration with my actual lab values. Knowing I'm APOE4 heterozygous would be curiosity. Knowing I'm APOE4 heterozygous with last quarter's LDL of whatever it was would be decision-relevant. The genome is one input. Lab results, current medications, and family history are the others. None of them mean as much in isolation as they do together. That's the next thing this tool needs to do.