From fcb854dca6fdaf59fe0753259d6cc8ed634b2650 Mon Sep 17 00:00:00 2001 From: JasonAHodgson Date: Sun, 4 Oct 2026 11:48:54 +0100 Subject: [PATCH 1/3] fix rbind() on loci with a missing allele_ref --- NEWS.md | 8 ++ R/rbind_dry_run.R | 8 +- tests/testthat/test_missing_alleles.R | 115 ++++++++++++++++++++++++++ 3 files changed, 130 insertions(+), 1 deletion(-) create mode 100644 tests/testthat/test_missing_alleles.R diff --git a/NEWS.md b/NEWS.md index 6616a56c..af1cf639 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,12 @@ # tidypopgen `dev` +* fix `rbind()` and `rbind_dry_run()` failing on a locus whose missing + allele is in `allele_ref`. NA alleles are replaced with "0" before + matching, but only `allele_alt` was covered; since `gen_tibble_bed()` maps + `allele_ref` to the bim file's allele2 column, a missing allele there left + an NA that reached a logical subscript and failed with "NAs are not + allowed in subscripted assignments" under `flip_strand = TRUE`. Such loci + are now dropped; a missing `allele_alt` is still resolved from the other + dataset as before. * implement an `autoplot()` method for distance matrices created by `pairwise_*` functions. This requires the matrices to have an appropriate class, so it will not work with objects generated by older versions of `tidypopgen` diff --git a/R/rbind_dry_run.R b/R/rbind_dry_run.R index 7f687468..2837481f 100644 --- a/R/rbind_dry_run.R +++ b/R/rbind_dry_run.R @@ -96,9 +96,15 @@ rbind_dry_run <- function( ref_df <- ref %>% show_loci() ref_df <- ref_df %>% mutate(id = seq_len(nrow(ref_df))) # replace NA with "0" for missing allele to avoid subsetting headaches - # (NA does not play nice with subsetting) + # (NA does not play nice with subsetting). Both allele columns need this: + # gen_tibble_bed() maps allele_ref to the bim file's allele2 column, so a + # missing allele can arrive in either one, and an NA left in allele_ref + # propagates through the == comparisons below into the logical vectors used + # as subscripts. ref_df$allele_alt[is.na(ref_df$allele_alt)] <- "0" target_df$allele_alt[is.na(target_df$allele_alt)] <- "0" + ref_df$allele_ref[is.na(ref_df$allele_ref)] <- "0" + target_df$allele_ref[is.na(target_df$allele_ref)] <- "0" # replace the names with a combination of chromosome and position if (use_position) { diff --git a/tests/testthat/test_missing_alleles.R b/tests/testthat/test_missing_alleles.R new file mode 100644 index 00000000..267575e1 --- /dev/null +++ b/tests/testthat/test_missing_alleles.R @@ -0,0 +1,115 @@ +# Loci with an unobserved allele. +# +# rbind_dry_run() replaces NA alleles with "0" before matching, because NA does +# not survive the logical subscripting used downstream. That substitution +# originally covered allele_alt only, so a missing allele arriving in allele_ref +# -- which is where gen_tibble_bed() puts the bim file's allele2 column -- +# left an NA in place and the call failed under flip_strand = TRUE with +# "NAs are not allowed in subscripted assignments". +# +# Note when extending this fixture: strand-ambiguous pairs (A/T, C/G) are +# removed when flip_strand = TRUE, and that applies to the alleles *after* any +# missing allele has been resolved from the other dataset. Every pair below is +# chosen to be non-ambiguous post-resolution. + +indiv_ref <- data.frame( + id = c("a", "b", "c"), population = c("pop1", "pop1", "pop2") +) +indiv_tgt <- data.frame( + id = c("x", "y", "z"), population = c("pop3", "pop3", "pop4") +) +geno <- rbind( + c(1, 1, 0, 1, 1, 0), + c(2, 1, 0, 0, 0, 0), + c(2, 2, 0, 0, 1, 1) +) +base_loci <- data.frame( + name = paste0("rs", 1:6), + chromosome = paste0("chr", c(1, 1, 1, 1, 2, 2)), + position = as.integer(c(3, 5, 65, 343, 23, 456)), + genetic_dist = as.double(rep(0, 6)) +) + +# ref target expected +# rs1 A/G A/G kept, matches as is +# rs2 T/C C/T kept, needs swap +# rs3 C/NA C/T kept, ref allele_alt resolved from target +# rs4 NA/G A/G dropped: gap is in allele_ref, not resolvable +# rs5 A/G G/NA kept, target allele_alt resolved, needs swap +# rs6 T/C A/G kept, needs strand flip +loci_ref <- cbind(base_loci, data.frame( + allele_ref = c("A", "T", "C", NA, "A", "T"), + allele_alt = c("G", "C", NA, "G", "G", "C") +)) +loci_tgt <- cbind(base_loci, data.frame( + allele_ref = c("A", "C", "C", "A", "G", "A"), + allele_alt = c("G", "T", "T", "G", NA, "G") +)) + +gt_ref <- gen_tibble( + x = geno, loci = loci_ref, indiv_meta = indiv_ref, + valid_alleles = c("A", "T", "C", "G"), quiet = TRUE +) +gt_tgt <- gen_tibble( + x = geno, loci = loci_tgt, indiv_meta = indiv_tgt, + valid_alleles = c("A", "T", "C", "G"), quiet = TRUE +) + +kept_loci <- c("rs1", "rs2", "rs3", "rs5", "rs6") + + +test_that("rbind_dry_run does not error on a missing allele in allele_ref", { + # the regression: before the fix, the NA left in allele_ref reached a logical + # subscript and the call failed + expect_no_error( + rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) + ) + expect_no_error( + rbind_dry_run(gt_ref, gt_tgt, flip_strand = FALSE, quiet = TRUE) + ) +}) + +test_that("no NA leaks into the merge report", { + report <- rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) + expect_false(anyNA(report$target$to_flip)) + expect_false(anyNA(report$target$to_swap)) + expect_false(anyNA(report$target$name)) + expect_false(anyNA(report$ref$name)) +}) + +test_that("a missing allele_alt is still resolved from the other dataset", { + report <- rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) + # rs3's gap is on the reference side, rs5's on the target side + expect_equal(report$ref$missing_allele[report$ref$name == "rs3"], "T") + expect_equal(report$target$missing_allele[report$target$name == "rs5"], "A") + expect_false(is.na(report$target$new_id[report$target$name == "rs3"])) + expect_false(is.na(report$target$new_id[report$target$name == "rs5"])) +}) + +test_that("a missing allele_ref is dropped rather than resolved", { + # resolve_missing_alleles() only inspects allele_1 (i.e. allele_alt), so a gap + # in allele_ref cannot be recovered. It must be dropped quietly, not error. + report <- rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) + expect_true(is.na(report$target$new_id[report$target$name == "rs4"])) + expect_false(report$target$to_flip[report$target$name == "rs4"]) + expect_false(report$target$to_swap[report$target$name == "rs4"]) +}) + +test_that("normal harmonisation is unaffected", { + report <- rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) + expect_true(report$target$to_swap[report$target$name == "rs2"]) + expect_true(report$target$to_swap[report$target$name == "rs5"]) + expect_true(report$target$to_flip[report$target$name == "rs6"]) + expect_setequal(report$target$name[!is.na(report$target$new_id)], kept_loci) +}) + +test_that("rbind produces a merged object with no missing alleles", { + merged <- rbind(gt_ref, gt_tgt, + flip_strand = TRUE, quiet = TRUE, + backingfile = tempfile("test_missing_allele_") + ) + expect_equal(nrow(merged), nrow(gt_ref) + nrow(gt_tgt)) + expect_setequal(loci_names(merged), kept_loci) + expect_false(anyNA(show_loci(merged)$allele_ref)) + expect_false(anyNA(show_loci(merged)$allele_alt)) +}) From 618ff51481b521f53e9b330cafe1ce739f4494c9 Mon Sep 17 00:00:00 2001 From: JasonAHodgson Date: Sun, 4 Oct 2026 12:51:58 +0100 Subject: [PATCH 2/3] fix rbind() on loci with a missing allele_ref --- tests/testthat/test_missing_alleles.R | 115 -------------------------- 1 file changed, 115 deletions(-) delete mode 100644 tests/testthat/test_missing_alleles.R diff --git a/tests/testthat/test_missing_alleles.R b/tests/testthat/test_missing_alleles.R deleted file mode 100644 index 267575e1..00000000 --- a/tests/testthat/test_missing_alleles.R +++ /dev/null @@ -1,115 +0,0 @@ -# Loci with an unobserved allele. -# -# rbind_dry_run() replaces NA alleles with "0" before matching, because NA does -# not survive the logical subscripting used downstream. That substitution -# originally covered allele_alt only, so a missing allele arriving in allele_ref -# -- which is where gen_tibble_bed() puts the bim file's allele2 column -- -# left an NA in place and the call failed under flip_strand = TRUE with -# "NAs are not allowed in subscripted assignments". -# -# Note when extending this fixture: strand-ambiguous pairs (A/T, C/G) are -# removed when flip_strand = TRUE, and that applies to the alleles *after* any -# missing allele has been resolved from the other dataset. Every pair below is -# chosen to be non-ambiguous post-resolution. - -indiv_ref <- data.frame( - id = c("a", "b", "c"), population = c("pop1", "pop1", "pop2") -) -indiv_tgt <- data.frame( - id = c("x", "y", "z"), population = c("pop3", "pop3", "pop4") -) -geno <- rbind( - c(1, 1, 0, 1, 1, 0), - c(2, 1, 0, 0, 0, 0), - c(2, 2, 0, 0, 1, 1) -) -base_loci <- data.frame( - name = paste0("rs", 1:6), - chromosome = paste0("chr", c(1, 1, 1, 1, 2, 2)), - position = as.integer(c(3, 5, 65, 343, 23, 456)), - genetic_dist = as.double(rep(0, 6)) -) - -# ref target expected -# rs1 A/G A/G kept, matches as is -# rs2 T/C C/T kept, needs swap -# rs3 C/NA C/T kept, ref allele_alt resolved from target -# rs4 NA/G A/G dropped: gap is in allele_ref, not resolvable -# rs5 A/G G/NA kept, target allele_alt resolved, needs swap -# rs6 T/C A/G kept, needs strand flip -loci_ref <- cbind(base_loci, data.frame( - allele_ref = c("A", "T", "C", NA, "A", "T"), - allele_alt = c("G", "C", NA, "G", "G", "C") -)) -loci_tgt <- cbind(base_loci, data.frame( - allele_ref = c("A", "C", "C", "A", "G", "A"), - allele_alt = c("G", "T", "T", "G", NA, "G") -)) - -gt_ref <- gen_tibble( - x = geno, loci = loci_ref, indiv_meta = indiv_ref, - valid_alleles = c("A", "T", "C", "G"), quiet = TRUE -) -gt_tgt <- gen_tibble( - x = geno, loci = loci_tgt, indiv_meta = indiv_tgt, - valid_alleles = c("A", "T", "C", "G"), quiet = TRUE -) - -kept_loci <- c("rs1", "rs2", "rs3", "rs5", "rs6") - - -test_that("rbind_dry_run does not error on a missing allele in allele_ref", { - # the regression: before the fix, the NA left in allele_ref reached a logical - # subscript and the call failed - expect_no_error( - rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) - ) - expect_no_error( - rbind_dry_run(gt_ref, gt_tgt, flip_strand = FALSE, quiet = TRUE) - ) -}) - -test_that("no NA leaks into the merge report", { - report <- rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) - expect_false(anyNA(report$target$to_flip)) - expect_false(anyNA(report$target$to_swap)) - expect_false(anyNA(report$target$name)) - expect_false(anyNA(report$ref$name)) -}) - -test_that("a missing allele_alt is still resolved from the other dataset", { - report <- rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) - # rs3's gap is on the reference side, rs5's on the target side - expect_equal(report$ref$missing_allele[report$ref$name == "rs3"], "T") - expect_equal(report$target$missing_allele[report$target$name == "rs5"], "A") - expect_false(is.na(report$target$new_id[report$target$name == "rs3"])) - expect_false(is.na(report$target$new_id[report$target$name == "rs5"])) -}) - -test_that("a missing allele_ref is dropped rather than resolved", { - # resolve_missing_alleles() only inspects allele_1 (i.e. allele_alt), so a gap - # in allele_ref cannot be recovered. It must be dropped quietly, not error. - report <- rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) - expect_true(is.na(report$target$new_id[report$target$name == "rs4"])) - expect_false(report$target$to_flip[report$target$name == "rs4"]) - expect_false(report$target$to_swap[report$target$name == "rs4"]) -}) - -test_that("normal harmonisation is unaffected", { - report <- rbind_dry_run(gt_ref, gt_tgt, flip_strand = TRUE, quiet = TRUE) - expect_true(report$target$to_swap[report$target$name == "rs2"]) - expect_true(report$target$to_swap[report$target$name == "rs5"]) - expect_true(report$target$to_flip[report$target$name == "rs6"]) - expect_setequal(report$target$name[!is.na(report$target$new_id)], kept_loci) -}) - -test_that("rbind produces a merged object with no missing alleles", { - merged <- rbind(gt_ref, gt_tgt, - flip_strand = TRUE, quiet = TRUE, - backingfile = tempfile("test_missing_allele_") - ) - expect_equal(nrow(merged), nrow(gt_ref) + nrow(gt_tgt)) - expect_setequal(loci_names(merged), kept_loci) - expect_false(anyNA(show_loci(merged)$allele_ref)) - expect_false(anyNA(show_loci(merged)$allele_alt)) -}) From 12e0f0e42e1c7a59ab2d14eddc1d8a1324d87902 Mon Sep 17 00:00:00 2001 From: JasonAHodgson Date: Sun, 4 Oct 2026 13:22:46 +0100 Subject: [PATCH 3/3] move missing-allele tests into test_rbind_dry_run.R --- tests/testthat/test_rbind_dry_run.R | 136 ++++++++++++++++++++++++++++ 1 file changed, 136 insertions(+) diff --git a/tests/testthat/test_rbind_dry_run.R b/tests/testthat/test_rbind_dry_run.R index 26f96a1e..aeee84fa 100644 --- a/tests/testthat/test_rbind_dry_run.R +++ b/tests/testthat/test_rbind_dry_run.R @@ -154,6 +154,142 @@ test_that("missing cases are given the correct alleles", { expect_true(miss_pop_a_flipped_swapped$to_flip) }) +# Loci with an unobserved allele ------------------------------------------ +# +# A locus can carry an unobserved allele on either side of a merge, and the +# gap can sit in either allele column: gen_tibble_bed() maps allele_ref to +# the bim file's allele2 column and allele_alt to allele1, so a missing code +# in either column of the source file can arrive in either column here. +# These loci must either be harmonised using the allele supplied by the +# other dataset, or dropped quietly -- never carried into the report as NA, +# which propagates into the logical subscripting used downstream. +# +# Note when extending this fixture: strand-ambiguous pairs (A/T, C/G) are +# removed when flip_strand = TRUE, and that applies to the alleles *after* +# any missing allele has been resolved from the other dataset. Every pair +# below is chosen to be non-ambiguous post-resolution. + +unobs_indiv_ref <- data.frame( + id = c("a", "b", "c"), population = c("pop1", "pop1", "pop2") +) +unobs_indiv_tgt <- data.frame( + id = c("x", "y", "z"), population = c("pop3", "pop3", "pop4") +) +unobs_geno <- rbind( + c(1, 1, 0, 1, 1, 0), + c(2, 1, 0, 0, 0, 0), + c(2, 2, 0, 0, 1, 1) +) +unobs_base_loci <- data.frame( + name = paste0("rs", 1:6), + chromosome = paste0("chr", c(1, 1, 1, 1, 2, 2)), + position = as.integer(c(3, 5, 65, 343, 23, 456)), + genetic_dist = as.double(rep(0, 6)) +) + +# ref target expected +# rs1 A/G A/G kept, matches as is +# rs2 T/C C/T kept, needs swap +# rs3 C/NA C/T kept, ref allele_alt resolved from target +# rs4 NA/G A/G dropped: gap is in allele_ref, not resolvable +# rs5 A/G G/NA kept, target allele_alt resolved, needs swap +# rs6 T/C A/G kept, needs strand flip +unobs_loci_ref <- cbind(unobs_base_loci, data.frame( + allele_ref = c("A", "T", "C", NA, "A", "T"), + allele_alt = c("G", "C", NA, "G", "G", "C") +)) +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") +)) + +unobs_gt_ref <- gen_tibble( + x = unobs_geno, loci = unobs_loci_ref, indiv_meta = unobs_indiv_ref, + valid_alleles = c("A", "T", "C", "G"), quiet = TRUE +) +unobs_gt_tgt <- gen_tibble( + x = unobs_geno, loci = unobs_loci_tgt, indiv_meta = unobs_indiv_tgt, + valid_alleles = c("A", "T", "C", "G"), quiet = TRUE +) + +unobs_kept_loci <- c("rs1", "rs2", "rs3", "rs5", "rs6") + + +test_that("an unobserved allele in either column does not error", { + # the gap sits in allele_ref for rs4 and in allele_alt for rs3 and rs5 + expect_no_error( + rbind_dry_run(unobs_gt_ref, unobs_gt_tgt, + flip_strand = TRUE, quiet = TRUE + ) + ) + expect_no_error( + rbind_dry_run(unobs_gt_ref, unobs_gt_tgt, + flip_strand = FALSE, quiet = TRUE + ) + ) +}) + +test_that("unobserved alleles leave no NA in the merge report", { + unobs_report <- rbind_dry_run(unobs_gt_ref, unobs_gt_tgt, + flip_strand = TRUE, quiet = TRUE + ) + expect_false(anyNA(unobs_report$target$to_flip)) + expect_false(anyNA(unobs_report$target$to_swap)) + expect_false(anyNA(unobs_report$target$name)) + expect_false(anyNA(unobs_report$ref$name)) +}) + +test_that("an unobserved allele_alt is resolved from the other dataset", { + unobs_report <- rbind_dry_run(unobs_gt_ref, unobs_gt_tgt, + flip_strand = TRUE, quiet = TRUE + ) + # rs3's gap is on the reference side, rs5's on the target side + tgt <- unobs_report$target + expect_equal(unobs_report$ref$missing_allele[ + unobs_report$ref$name == "rs3" + ], "T") + expect_equal(tgt$missing_allele[tgt$name == "rs5"], "A") + expect_false(is.na(tgt$new_id[tgt$name == "rs3"])) + expect_false(is.na(tgt$new_id[tgt$name == "rs5"])) +}) + +test_that("an unobserved allele_ref is dropped rather than resolved", { + # resolve_missing_alleles() only inspects allele_1 (i.e. allele_alt), so a + # gap in allele_ref cannot be recovered from the other dataset. Such a + # locus is dropped from the merge without error. + unobs_report <- rbind_dry_run(unobs_gt_ref, unobs_gt_tgt, + flip_strand = TRUE, quiet = TRUE + ) + tgt <- unobs_report$target + expect_true(is.na(tgt$new_id[tgt$name == "rs4"])) + expect_false(tgt$to_flip[tgt$name == "rs4"]) + expect_false(tgt$to_swap[tgt$name == "rs4"]) +}) + +test_that("swaps and flips are still detected alongside unobserved alleles", { + unobs_report <- rbind_dry_run(unobs_gt_ref, unobs_gt_tgt, + flip_strand = TRUE, quiet = TRUE + ) + tgt <- unobs_report$target + expect_true(tgt$to_swap[tgt$name == "rs2"]) + expect_true(tgt$to_swap[tgt$name == "rs5"]) + expect_true(tgt$to_flip[tgt$name == "rs6"]) + expect_setequal(tgt$name[!is.na(tgt$new_id)], unobs_kept_loci) +}) + +test_that("rbind of loci with unobserved alleles gives a complete object", { + merged <- rbind(unobs_gt_ref, unobs_gt_tgt, + flip_strand = TRUE, quiet = TRUE, + backingfile = tempfile("test_unobserved_allele_") + ) + expect_equal(nrow(merged), nrow(unobs_gt_ref) + nrow(unobs_gt_tgt)) + expect_setequal(loci_names(merged), unobs_kept_loci) + expect_false(anyNA(show_loci(merged)$allele_ref)) + expect_false(anyNA(show_loci(merged)$allele_alt)) +}) + + + # # # #reference file reordered