Skip to content

fix rbind() on loci with a missing allele_ref - #207

Open
JasonAHodgson wants to merge 3 commits into
devfrom
fix-missing-allele-rbind
Open

JasonAHodgson wants to merge 3 commits into
devfrom
fix-missing-allele-rbind

Conversation

@JasonAHodgson

Copy link
Copy Markdown

Fix rbind() erroring on a locus whose missing allele is in allele_ref

The bug

rbind_dry_run(hgdp_gt, tg_gt, flip_strand = TRUE, quiet = TRUE)
#> Error in target_sub$allele_1[!to_keep_orig] <- flip(target_sub$allele_1[!to_keep_orig]) :
#>   NAs are not allowed in subscripted assignments

rbind_dry_run() already guards against this, at R/rbind_dry_run.R:99:

# replace NA with "0" for missing allele to avoid subsetting headaches
# (NA does not play nice with subsetting)
ref_df$allele_alt[is.na(ref_df$allele_alt)] <- "0"
target_df$allele_alt[is.na(target_df$allele_alt)] <- "0"

but only for allele_alt. gen_tibble_bed() maps

allele_ref = as.character(bim$allele2),
allele_alt = as.character(bim$allele1)

so a missing allele in the .bim file's allele2 column becomes a missing
allele_ref, keeps its NA through harmonise_missing_values(), and reaches the
== comparisons that build to_keep_orig. NA propagates, !to_keep_orig is
NA, and the subscripted assignment fails.

The fix

Sanitise both allele columns:

ref_df$allele_ref[is.na(ref_df$allele_ref)] <- "0"
target_df$allele_ref[is.na(target_df$allele_ref)] <- "0"

Nothing else changes. Loci with a missing allele_alt are still recovered by
resolve_missing_alleles() as before; loci whose gap is in allele_ref are now
dropped quietly instead of erroring.

How it came up

Merging HGDP genotyped on the Affymetrix Axiom Human Origins array with 1000
Genomes phase 3. The Axiom .bim writes X in the allele2 column for 96 loci
with no observed second allele. Declaring it via the documented route —

gen_tibble("HGDP_Axiom_allmap.bed", missing_alleles = c("0", ".", "X"), ...)

— produces an object rbind() cannot merge.

Ordinary plink --make-bed output is mostly unaffected: PLINK writes 0 in the
A1 column, which maps to allele_alt and was always sanitised. So this needs a
missing allele specifically in allele2, which is unusual but not exotic.

Tests

tests/testthat/test_missing_alleles.R, six tests over a six-locus fixture:

locus ref target covers
rs1 A/G A/G matches unchanged
rs2 T/C C/T swap still detected
rs3 C/NA C/T allele_alt gap on the ref side, resolved
rs4 NA/G A/G the regression — allele_ref gap, dropped not errored
rs5 A/G G/NA allele_alt gap on the target side, resolved
rs6 T/C A/G strand flip still detected

A warning for anyone extending this fixture, which is in a comment at the top of
the file: strand-ambiguous pairs (A/T, C/G) are dropped when
flip_strand = TRUE, and that applies to the alleles after a missing one has
been resolved from the other dataset. A pair that looks fine as written can
become ambiguous once repaired, and the locus then disappears for a reason that
has nothing to do with missing data.

Two things noticed but not changed here

  1. resolve_missing_alleles() only inspects allele_1 (i.e. allele_alt),
    so a gap in allele_ref is dropped even where the other dataset could supply
    the allele. Making it consider both columns would recover those loci. Left
    out to keep this diff to the crash.

  2. flip() is NA-unsafe. bases == "T" is NA for NA input, so
    bases[NA] <- "A" errors. Unreachable now that no NA survives to that
    point, but worth hardening — a lookup (c(A = "T", T = "A", C = "G", G = "C"))
    would be NA-safe and shorter, at the cost of needing care to preserve the
    documented pass-through of non-ACGT values like "m".

Happy to do either in a follow-up.

Closes #

@codecov

codecov Bot commented Oct 4, 2026 •

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 93.34%. Comparing base (406f87f) to head (12e0f0e).

Additional details and impacted files
@@           Coverage Diff           @@
##              dev     #207   +/-   ##
=======================================
  Coverage   93.34%   93.34%           
=======================================
  Files         133      133           
  Lines        7920     7922    +2     
=======================================
+ Hits         7393     7395    +2     
  Misses        527      527           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@dramanica

Copy link
Copy Markdown
Member

This looks good. However, can you please move the tests in test_rbind_dry_run.R? Also, make sure that the description and comments in tests do not refer to the bug, but rather are generic about what the tests check (so that, we someone looks at them in the future, they make sense). Thanks!

Copilot AI left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copilot review overview

🟡 Changes recommended

Shared sentinels can incorrectly retain loci with unresolved gaps, leaving missing allele metadata.

Review effort: Balanced
Findings: 1 High severity · 1 Medium severity

Open (2)
What changed in this PR

Fixes merge failures when allele_ref is missing by normalizing missing alleles before comparison.

Changes:

  • Normalizes missing reference alleles.
  • Adds regression coverage for missing alleles.
  • Documents the fix.
File Description
R/​rbind_dry_run.R Normalizes missing reference alleles.
tests/​testthat/​test_rbind_dry_run.R Adds merge regression tests.
NEWS.md Documents corrected behavior.

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

Comment thread R/rbind_dry_run.R
Comment on lines +106 to +107
ref_df$allele_ref[is.na(ref_df$allele_ref)] <- "0"
target_df$allele_ref[is.na(target_df$allele_ref)] <- "0"
Comment on lines +201 to +203
unobs_loci_tgt <- cbind(unobs_base_loci, data.frame(
allele_ref = c("A", "C", "C", "A", "G", "A"),
allele_alt = c("G", "T", "T", "G", NA, "G")

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants