Skip to content

Fix filtering bugs2 - #437

Merged
jeromekelleher merged 8 commits into
sgkit-dev:mainfrom
jeromekelleher:fix-filtering-bugs2
Jun 8, 2026
Merged

Fix filtering bugs2#437
jeromekelleher merged 8 commits into
sgkit-dev:mainfrom
jeromekelleher:fix-filtering-bugs2

Conversation

@jeromekelleher

Copy link
Copy Markdown
Member

Closes #435
Closes #436

The variant-scope logical operators & and | used bare np.logical_and/or,
which cannot align a 1-D per-variant mask (e.g. N_ALT) with a 2-D
per-allele or per-sample mask (e.g. INFO/AC, FMT/DP), raising a broadcast
ValueError. The && and || operators already handled this via dimension
expansion. Factor that expansion into _align_mask_dims and route & and |
through single_and/single_or so a 1-D operand broadcasts across the
allele/sample axis before the root collapse, matching bcftools.

This lets the PLINK conversion example use a single N_ALT == 1 & AC > 0
filter with one materialise pass instead of two.
Filter expressions comparing a per-allele field (Number=A, e.g. INFO/AC,
shape (variants, alt_alleles)) against a per-variant field (e.g. INFO/AN,
shape (variants,)) raised "operands could not be broadcast together"
because ComparisonOperator.eval and the BinaryOperator arithmetic path
fed both operands to numpy without aligning dimensions. Only the logical
operators aligned a 1-D operand to (n, 1). This blocked the canonical
bcftools idiom for keeping polymorphic sites, AC>0 && AC<AN.

Reuse the (renamed) _align_dims helper in the comparison and arithmetic
paths so a 1-D operand broadcasts across the allele axis.

Once tag-vs-tag comparisons evaluate, the missing-data handling also has
to match bcftools, which treats two missing values as equal: == is true
and != is false when both operands are absent, a present value never
equals a missing one, and ordering comparisons are false on any missing
operand. The previous single-fill model returned the opposite for the
both-missing case; it was only correct for tag-vs-constant comparisons.

Update the PLINK export docs to use AC>0 && AC<AN, which (unlike AC>0
alone) also drops sites monomorphic for ALT where AC==AN.
For the per-allele INFO/AC field, `&` requires the same allele to satisfy
both AC>0 and AC<AN, whereas `&&` would pass a multi-allelic site where the
two conditions are met by different alleles. Match format_conversion.md,
which already uses `&`.
A logical operator combining a per-allele condition (variant-scope 2-D,
axis 1 = alleles) with a per-sample condition (sample-scope 2-D, axis 1 =
samples) crashed with a numpy broadcast error, e.g. 'AC>0 & FMT/DP>0' — a
filter bcftools evaluates. BinaryOperator now collapses the allele axis of
each variant-scope operand to a per-variant decision before combining with a
per-sample mask, gated on the whole node being sample-scope so all-variant
expressions keep their same-allele semantics.

A per-allele vs per-sample comparison (INFO/AC > FMT/DP) is ill-defined and
now raises a clear error instead of a raw numpy ValueError. Scalar operands to
&& / || broadcast consistently with & / | rather than crashing.

Adds direct unit tests for the combination helpers (_align_dims, single_and,
double_and, _collapse_to_variant, ...), a dimension/scope/ragged-field matrix,
bcftools-validated mixed INFO/FORMAT cases, and the previously-uncovered
ParseError / no-op-filter / special-operator-repr paths, bringing
bcftools_filter.py to 100% line and branch coverage.
The filter pipeline used a bespoke _missing_mask that lumped INT_FILL
(Number=A trailing padding for variants with fewer alleles than the chunk
maximum) together with genuine missing values. Because a variant-scoped
per-allele result is collapsed at the root with np.any over the allele
axis, a fill slot forced to True leaked through and wrongly included the
variant. Verified against bcftools: `AC!=3` on a biallelic site whose only
allele is 3 was wrongly kept (the FILL pad made != True).

Replace _missing_mask with the precise utils.is_missing / utils.is_fill
helpers, splitting the node-tree predicate into is_missing() and is_fill()
(both propagated through arithmetic). ComparisonOperator drives bcftools
missing semantics from is_missing, then forces fill slots to False for
every operator so a non-existent allele never matches. FLAG fields carry
no missing mask (DB=0 must still match unset flags).
Add prose to the PLINK conversion example explaining that
materialise_variant_filter makes a full pass over the variant axis,
and enable INFO logging in the example so the reader can see and
interpret the passes over the data.
The dimension logic assumed a variant-scope 2-D operand was always
alt-allele-dimensioned (Number=A). INFO fields are 2-D for other reasons
(Number=R, fixed-N, Number=.) and those already collapsed correctly, but
combining two variant-scope 2-D fields of *different* widths under a logical
operator (e.g. 'IIR>0 | IID>0', which bcftools evaluates) crashed with a numpy
broadcast error.

_prepare_logical_operands now collapses each variant-scope 2-D operand to a
per-variant decision (any element) whenever the operands don't share one
element axis -- mixed scope, or differing per-allele widths -- while same-width
same-scope operands keep the element axis (so 'AC>0 & AC<2' stays a same-allele
test). Comparisons and arithmetic between genuinely incompatible-width 2-D
operands now raise a clear error instead of a raw numpy crash.

Tests exercise the full INFO/FORMAT Number matrix via field_type_combos.vcf.gz:
multi-valued INFO (Number=2/R/.) comparisons validated against bcftools,
different-width logical combos validated against bcftools, and multi-valued
(3-D) FORMAT fields confirmed to raise UnsupportedHigherDimensionalFormatFieldsError
across every cardinality and dtype. bcftools_filter.py stays at 100% coverage.
@jeromekelleher
jeromekelleher merged commit f99fde3 into sgkit-dev:main Jun 8, 2026
13 checks passed
@jeromekelleher
jeromekelleher deleted the fix-filtering-bugs2 branch June 8, 2026 10:37
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

1 participant