A1/A2 conventions for plink output #327
jeromekelleher
started this conversation in
General
Replies: 0 comments
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Uh oh!
There was an error while loading. Please reload this page.
Compatibility report — A1/A2 convention for plink-1 output
Question
vcztools currently writes
.bed/.bimusing plink 1.9's minor-allele convention (A2 = REF, A1 = the most-frequent non-REF allele) computed byWriter._compute_allelesinvcztools/plink.py:106–143. Computing this requires a full pass overcall_genotypeper variant.The plink-2-style REF/ALT ordering (A1 = ALT, A2 = REF) as the desired default. This report addresses the deeper question: can we drop the minor-allele path entirely, or must we keep it as an opt-in for some consumers?
TL;DR
Yes — minor-allele can be dropped entirely. None of the consumers in the IO study (admixture, bolt, flashpca, gcta, king, regenie, plink 1.9, plink 2) require A1 to be the minor allele as a semantic precondition. The convention is plink 1.9's historical default, not a correctness assumption. Three caveats — see "Risks" below.
Per-consumer findings
For each tool, the question is: does the tool produce wrong output if A1 is the ALT allele rather than the minor allele?
plink 1.9
--keep-allele-order(or--real-ref-alleles). The on-disk bytes never change — only the in-memory association.--freq,--assoc,--logistic, etc.) reflect the in-memory ordering, so an unflagged plink 1.9 run on our output gets minor-allele A1 in its outputs regardless of what we wrote on disk.--keep-allele-ordersee exactly the labelling we wrote; pipelines that don't get plink-1.9's frequency-reordered view. Both are correct interpretations of the same.bedbytes.plink 2
.bimformat spec (cog-genomics) literally defines A1 = ALT, A2 = REF.REGENIE
.bedinput: A2 is treated as the reference allele (so A1 is the effect/ALT). Confirmed by REGENIE docs:--ref-firstexists to opt out into "first allele = reference".BOLT-LMM
.bim.GCTA
--make-grm,--reml,--pca,--make-grm-bin, etc.) standardise genotypes per-variant before use. Standardisation is invariant to allele labelling (the variance is the same whether you count copies of A1 or A2).--cojoand friends consume summary statistics with explicit allele columns — orthogonal to.bedencoding.)KING
flashpca
--project) requires consistent allele coding between training and projection sets — but any consistent coding works; REF/ALT is in fact more stable across cohorts than minor-allele (the latter can flip between cohorts if a variant sits near 50% MAF).ADMIXTURE
Cross-cutting observations
--keep-allele-order. So even the "stay compatible with plink 1.9" framing dissolves: under default plink 1.9 invocations, the minor-allele relabelling happens at plink load time regardless of what we wrote.Argument for dropping minor-allele entirely
bincountover2 × num_sampleshaploid calls. For the streaming source this collapses init time from O(N × S) to O(N); forwrite_plinkit's a smaller win (genotypes are read for encoding anyway) but still a non-trivial reduction in arithmetic._a12_filled, the per-chunk write-through cache,eager_a12parameter ambiguity, and the §7 "internal design notes" complexity in the streaming spec. The a12 array becomes a static(N, 2) int8derived fromvariant_alleleat construction.compute_a12_minor_alleleas an opt-in means: two code paths, two test classes (TestComputeAllelescovers the algorithm), divergent benchmarks, mode plumbing throughwrite_plink, the streaming source, andview-plink1's eventual flag set.Argument for keeping minor-allele as opt-in
.bedwhose bytes encode minor-allele homozygotes as00. Re-running our pipeline on the same source VCF should ideally produce the same bytes. Counter: this kind of byte-level reproduction is brittle anyway (depends on plink 1.9's tie-breaking, dataset-specific frequencies, and same-cohort assumption). A semantically equivalent output is achievable viaplink1.9 --bfile our_out --keep-allele-order --recode bedfollowed by re-running plink 1.9 without that flag, which yields the legacy bytes.vcztoolsusers runningwrite_plinkmay have downstream comparisons assuming the current bytes. Counter: document the change in release notes and bump a minor version.Risks if we drop it
.bedbytes to a stored plink-1.9-generated.bed(rather than to genotype values) will see a diff. The fix is documented and one flag away..beddirectly and divide by 2 to get "MAF" without consulting.bim) could silently misinterpret. The.bimallele columns are still correct, so any tool that consults them gets the right answer. This is a script-bug class, not a vcztools-bug class.Recommendation
Drop minor-allele. The compatibility surface is empty for the consumers we care about, the performance and design wins are substantial, and the residual risks are document-able. Add a clear release-notes entry pointing users at
--keep-allele-orderfor plink 1.9 reading and explaining that REF/ALT is the on-disk convention.If the maintainer wants a softer landing, the alternative is a single release where
write_plink(..., a12_strategy="minor_allele")is accepted but emits aDeprecationWarning, with removal one release later. This is more friction for ~zero ecosystem benefit given the analysis above, so I'd skip it.Sources
All reactions