Part 1 of 5.
Thank you for reading this post, don't forget to subscribe!Every genotyping array tells a small story before it tells you anything about DNA. A glass slide goes into a scanner, the scanner writes two IDAT files per sample (one red, one green), and a few hundred thousand bead intensities wait for someone to turn them into genotypes. That someone is almost always a piece of closed software. You click, it thinks, and out come the calls.
For most people this is perfectly fine. The calls are good. The labs are busy. Nobody has time to ask how a pair of numbers became “AB”.
I had the time, and I had the question. SimurgArray is the answer I ended up writing.
The black box in the middle
Illumina Infinium arrays are everywhere. Biobanks use them, consumer genetics companies use them, research groups use them because they are cheap per sample and wonderfully boring in the best sense: they work. Hundreds of millions of people have been genotyped this way.
The pipeline is conceptually simple:
- Read the raw intensities from the IDAT files.
- Normalize them, because the red and green channels never quite agree with each other.
- Compare each sample to a cluster file (an EGT) that says where AA, AB and BB usually sit.
- Call genotypes, give each call a score, and write everything to a GTC file.
- Turn that into something the rest of the world reads, usually a VCF.
Steps 1 to 4 are where the interesting decisions live, and they have historically been done by software you can run but cannot read. Good open tools exist around the edges, and I will talk about them in the next post, but I could not find one place where the whole chain was open, tested against the vendor’s own outputs, and documented well enough that a stranger could check every step.
So I built that place.
What SimurgArray actually is
SimurgArray is a personal research toolkit that takes IDAT files all the way to QC reports, VCF and PLINK files, and then keeps going:
- Genotype calling that reproduces the vendor’s GTC files byte for byte in 2138 of 2141 public 1000 Genomes Omni2.5 samples. The last three are documented, not hidden.
- Sample and control QC with an HTML report you can actually click around in.
- Copy number, loss of heterozygosity and mosaicism calling across the genome.
- Pharmacogenomics: targeted copy number for genes like CYP2D6, star allele diplotypes and CPIC phenotypes.
- Methylation array control QC.
- ArrayTrain, which builds your own cluster file from your own cohort, so you are not forever borrowing someone else’s.
Everything is specified item by item in a functional spec with 416 numbered requirements, and every algorithm points to where it came from: a paper, a public technical note, an expired patent, or an open source project with a compatible license. When our numbers differ from a reference tool, the difference goes into a file called known_unknowns.md with a date and an explanation. Thirty entries so far. Some of them are my favorite reading.
Reason one: you should be able to check
Science has a slightly awkward relationship with tools it cannot inspect. We publish methods sections that say “genotypes were called with the manufacturer’s software, default settings”, and everybody nods, because what else would you say?
The problem shows up later. A cluster file that does not match your chips. A normalization version that changed between releases. A sex check that disagrees with the pedigree in 12 samples out of two thousand. With a black box you can only shrug. With open code you can open the file, find the line, and write down what happened.
SimurgArray was built around that idea. It has a differential testing harness that compares our outputs with reference tools field by field and labels each difference as exact, near, discrepancy or unknown. A nightly job runs it again at 02:30, when the Raspberry Pi has nothing better to do and the rest of the house is asleep.
Reason two: the cluster file problem
Here is something that surprised me more than it should have. Genotype calls depend heavily on the cluster file, and the cluster file was trained by someone else, on samples you have never seen, possibly for a slightly different version of your chip.
On a public dataset of 360 GSA arrays I worked with, the only public cluster file available was made for a sibling chip. It works, mostly. The quality metrics politely complain the whole time.
With ArrayTrain you can train a cluster file from your own cohort. On that dataset, a cluster file trained from 120 arrays gave 99.44 % concordance with 30x whole genome sequencing, against 99.46 % for the vendor’s file, while calling more genotypes. That is not a victory lap. It is a proof that the dependency is optional, which is a different and, I think, more useful thing.
Reason three: it should run on a small computer
A lot of genomics software quietly assumes a server room. SimurgArray was developed and mostly validated on a Raspberry Pi 5 with 8 GB of memory. Heavy jobs run inside a memory guard, because I learned the hard way what happens when the kernel has to choose between my test suite and my web browser. The browser lost. So did the test suite, actually.
With the optional Rust kernels a GSA sample takes about a second on a single Pi core. On a cloud server core it is roughly 80 samples a minute. The few cloud validation runs I needed cost a few dollars in total, and every one of them was set to switch itself off and was checked afterwards to make sure nothing kept billing quietly in the background.
I care about this because the people who would benefit most from open array analysis are often the ones with the least infrastructure. More on that in part 3.
Reason four: understanding is a feature
There is a kind of knowledge you only get by rebuilding something. I now know why the normalization uses the exact float operations it uses, why a median of an even number of values is not always the median you expect, and why a NaN can have a sign depending on which processor produced it. None of this is glamorous. All of it matters if you want two programs to agree to the last bit.
That understanding is now written down, in code and in documents, where it can outlive my own memory of it. Given how my memory works, that is a real advantage.
What it is not
SimurgArray is a research tool. It is not a medical device, it is not for diagnosis, and its pharmacogenomics output is not a prescription. I say this plainly because PGx reports look very official, and a well formatted table can be more convincing than it deserves to be.
It also does not use anything it should not. No decompiled binaries, no vendor data in the repository, no copied GPL code. Tools under the GPL are only ever called as separate programs to compare against. When I borrowed from MIT, BSD or Apache licensed projects, the credit is in THIRD_PARTY_NOTICES.md, word for word.
The short version
SimurgArray exists because array genotyping is too widely used to be understood by only one company, because cluster files should be something you can make yourself, and because a curious person with a small computer should be able to follow every number from the scanner to the VCF.
In the next post I will look at what else is out there, which is quite a lot, and why I still decided to write my own. The honest answer involves both good reasons and a certain amount of stubbornness.