Error in `downsample_immundata(idata, n = 2, seed = 33)`: No barcode units were selected. Increase `n`.

Expected `all(ann_na_counts == 0)` to be TRUE.
Differences:
`actual`:   FALSE
`expected`: TRUE 

context: single-chain duplicated receptors
annotation shape: 5x16
repertoire shape: 2x4
missing annotation columns: <none>
missing repertoire columns: <none>
annotation NA counts: imd_receptor_id=0, imd_repertoire_id=5, imd_count=5, imd_proportion=5, n_repertoires=0
repertoire NA counts: imd_repertoire_id=0, n_barcodes=0, n_receptors=2
unmatched repertoire ids: NA

annotations dump:
# A tibble: 5 x 16
  cell_id sample_id v_call j_call junction_aa locus umi_count filename          
* <chr>   <chr>     <chr>  <chr>  <chr>       <chr>     <dbl> <chr>             
1 cell1   Sample1   IGHV1  IGHJ1  CARW        IGH         100 "D:\\temp\\2026_0~
2 cell2   Sample1   IGHV1  IGHJ1  CARW        IGH         150 "D:\\temp\\2026_0~
3 cell3   Sample1   IGHV2  IGHJ2  CBRW        IGH         200 "D:\\temp\\2026_0~
4 cell4   Sample2   IGHV3  IGHJ3  CCRW        IGH         250 "D:\\temp\\2026_0~
5 cell5   Sample2   IGHV4  IGHJ4  CDRW        IGH         300 "D:\\temp\\2026_0~
# i 8 more variables: imd_barcode <chr>, imd_chain_id <int>,
#   imd_receptor_id <int>, imd_n_chains <dbl>, imd_count <dbl>,
#   imd_repertoire_id <int>, imd_proportion <dbl>, n_repertoires <int>

repertoires dump:
# A tibble: 2 x 4
  imd_repertoire_id sample_id n_barcodes n_receptors
              <int> <chr>          <dbl>       <int>
1                 1 Sample2            2          NA
2                 2 Sample1            3          NA

# TEST 1: Single-chain data with duplicated receptors
test_that("agg_repertoires counts single-chain receptors correctly", {
# Sample1: 3 cells, 2 unique receptors (cells 1&2 share receptor)
cell_id = c("cell1", "cell2", "cell3", "cell4", "cell5"),
sample_id = c("Sample1", "Sample1", "Sample1", "Sample2", "Sample2"),
v_call = c("IGHV1", "IGHV1", "IGHV2", "IGHV3", "IGHV4"),
j_call = c("IGHJ1", "IGHJ1", "IGHJ2", "IGHJ3", "IGHJ4"),
junction_aa = c("CARW", "CARW", "CBRW", "CCRW", "CDRW"), # cell1 & cell2 identical
locus = c("IGH", "IGH", "IGH", "IGH", "IGH"),
idata_agg <- agg_repertoires(idata, schema = "sample_id")
repertoires <- idata_agg$repertoires |> collect()
annotations <- idata_agg$annotations |> collect()
# Sample1: n_barcodes = 3, n_receptors = 2 (cells 1&2 share receptor)
# Sample2: n_barcodes = 2, n_receptors = 2
sample1_stats <- repertoires |> filter(sample_id == "Sample1")
sample2_stats <- repertoires |> filter(sample_id == "Sample2")
info = "Sample1 should have 2 unique receptors"
info = "Sample2 should have 2 unique receptors"
# Verify receptor counts within repertoires
# The shared receptor in Sample1 should have count = 2
info = "Shared receptor should have count of 2"
# TEST 2: Paired-chain data - potential double counting issue
test_that("agg_repertoires identifies double-counting in paired-chain data", {
# Each receptor appears TWICE (once per chain)
cell_id = c("cell1", "cell1", "cell2", "cell2", "cell3", "cell3"),
sample_id = c("Sample1", "Sample1", "Sample1", "Sample1", "Sample2", "Sample2"),
v_call = c("IGHV1", "IGLV1", "IGHV2", "IGLV2", "IGHV3", "IGLV3"),
j_call = c("IGHJ1", "IGLJ1", "IGHJ2", "IGLJ2", "IGHJ3", "IGLJ3"),
junction_aa = c("CARW", "CASL", "CBRW", "CBSL", "CCRW", "CCSL"),
locus = c("IGH", "IGL", "IGH", "IGL", "IGH", "IGL"),
umi_count = c(100, 100, 150, 150, 200, 200) # Same count for both chains
# Current behavior (with double-counting bug):
# Sample1: n_barcodes will be 4 (2 cells × 2 chains) - WRONG!
# Sample2: n_barcodes will be 2 (1 cell × 2 chains) - WRONG!
# Sample1: n_barcodes = 2 (actual cells)
# These tests document the CURRENT behavior (with double-counting)
# They should FAIL when the bug is fixed
# Test for n_receptors (this should be correct)
info = "Sample2 should have 1 unique receptor"
# TEST 3: Mixed scenario - paired-chain with shared receptors
test_that("agg_repertoires handles paired-chain data with shared receptors", {
# Create paired-chain data where some cells share receptors
sample_id = c("Sample1", "Sample1", "Sample1", "Sample1", "Sample1", "Sample1"),
v_call = c("IGHV1", "IGLV1", "IGHV1", "IGLV1", "IGHV2", "IGLV2"),
j_call = c("IGHJ1", "IGLJ1", "IGHJ1", "IGLJ1", "IGHJ2", "IGLJ2"),
junction_aa = c("CARW", "CASL", "CARW", "CASL", "CBRW", "CBSL"),
# - 2 unique receptors (cells 1&2 share)
# Current (incorrect due to double-counting):
# - n_barcodes will show 6 (3 cells × 2 chains)
# The shared receptor should have count = 2 × 2 = 4 (double-counted)
# Should be 2 (two cells sharing the receptor)
# TEST 4: Proportions calculation with double-counting
test_that("agg_repertoires proportions are affected by double-counting", {
# Simple paired-chain data for clear proportion calculation
cell_id = c("cell1", "cell1", "cell2", "cell2"),
sample_id = c("Sample1", "Sample1", "Sample1", "Sample1"),
v_call = c("IGHV1", "IGLV1", "IGHV2", "IGLV2"),
j_call = c("IGHJ1", "IGLJ1", "IGHJ2", "IGLJ2"),
junction_aa = c("CARW", "CASL", "CBRW", "CBSL"),
umi_count = c(100, 100, 100, 100) # Equal counts
# TEST 5: Verify n_repertoires calculation (should be correct)
test_that("agg_repertoires n_repertoires calculation is correct for paired data", {
# Create data where one receptor appears in multiple samples
sample_id = c("Sample1", "Sample1", "Sample2", "Sample2", "Sample3", "Sample3"),
# Check n_repertoires for the shared receptor
# First receptor appears in Sample1 and Sample2 (n_repertoires = 2)
# Second receptor appears only in Sample3 (n_repertoires = 1)
info = "Shared receptor should appear in 2 repertoires"
info = "Unique receptor should appear in 1 repertoire"
test_that("agg_repertoires preserves second chain data (no NAs)", {
# The input had 2 rows (1 cell x 2 chains).
# (The bug caused this to drop to 1 row).
info = "Annotation table should retain rows for both chains (IGH and IGL)"
# Filter for the Light Chain row and ensure V-call is present (not NA)
igl_row <- annotations |> filter(locus == "IGL")
info = "Should find exactly one row for the IGL chain"
info = "IGL row should have valid v_call data (not NA)"
info = "IGL v_call should match input data"
# Ensure the calculated stats (which are per-receptor) are mapped
# Both IGH and IGL rows belong to the same receptor, so they should
# both have the same imd_count and imd_proportion.
igh_row <- annotations |> filter(locus == "IGH")
info = "Both chains of the same receptor should share the same count stats"
info = "Both chains should belong to the same repertoire ID"