mpileup: use 64-bit accumulators in the Mann-Whitney bias Z scores

[cindykrafft]

Fixes #2594.

calc_mwu_biasZ() summed the tie adjustment (p*p-1)*p, the pair counts and the
group sizes in int. The tie term is cubic in the number of reads sharing one quality
or position bin and exceeds INT_MAX once a bin holds 1291 reads, which happens at a
single deep site (-d raised for amplicon, mitochondrial or viral data) and, under the
default -d 250, across the samples of a multi-sample call, since the histograms are
pooled over all samples. The wrapped term inflated the variance and shrank
INFO/MQBZ, BQBZ, RPBZ, SCBZ, MQSBZ and NMBZ towards zero: 1300 MAPQ-60
reference reads and 40 MAPQ-30 alternate reads give MQBZ=-7.88 where the
tie-corrected Z is −36.59 (1290 reference reads give the right −36.46).

This changes the accumulators to int64_t (no change in output where nothing
overflowed: make test otherwise unchanged) and adds a regression test,
test/mpileup/mwu-biasZ.1.bam (1,340 20-bp reads, 4.9 kB) with its reference and
expected output, run through test/test.pl as -B -d 100000 -a -AD -r ref:150; the
test fails on unmodified develop (MQBZ=-7.88325) and passes with the change.

NEWS entry included under the unreleased heading.

Assisted-by: Claude:claude-fable-5-1

Found in Mytochondria, a volunteer project that checks the numerical core of research software and verifies every finding by execution (methods and harnesses: https://github.com/cindykrafft/mytochondria/tree/main/audits/bcftools)

---

Generated by Claude Code

[jkbonfield]

Could you paste in the mcve_mpileup_mqbz_overflow.sh script referred to in the issue please? I tried reproducing it and I can get overflows, but much later than 1291. I'll experiment some more.

I'm not sure MQ is the first to overflow either, but it's an easy one to target as it's a single value. (In theory MQUAL can go up to 254, although in practice it's normally capped.)

The overflow comes from the cube root of 2^31. For 2^63 that's about 2.1 million. Hence still theoretically possible with 64-bit, but extreme and not something to be overly concerned with I feel. However I also asked Claude if it could rearrange the formulae to avoid accumulating cubed values which is promptly did and gave me an alternative that may also solve some other issues:

"Here's a refactor that removes the overflow and fixes a few related ones in the same function. The current develop code also overflows e, l, nanb and (na+nb)(na+nb-1) as int. The last of these is reachable in big multi-sample runs once N passes about 46k. This is the bug reported in #2594.

The idea: don't compute sum(p³) at all. The tie-corrected variance is

var = na·nb/12 · ((N+1) − Σ(p³−p) / (N(N−1)))

Because Σp = N, this simplifies to

var = na·nb · (N³ − Σp³) / (12·N·(N−1))

You can build N³ − Σp³ one bin at a time using only non-negative terms. Adding a bin of size p to a running total S contributes (S+p)³ − S³ − p³ = 3·S·p·(S+p). That gives three benefits:

Nothing gets cubed in integer arithmetic, so there's no overflow.
It avoids the cancellation in (N+1) - t/(N(N-1)), which subtracts two nearly equal large numbers.
The "all values in one bin" case gives exactly 0, not a tiny positive rounding error that would blow up Z."

I'm still evaluating whether this is a better alternative. For reference, the code it gave me was this:

double calc_mwu_biasZ(int *a, int *b, int n, int left_only, int do_Z) {
    int i;

    // Optimisation: with no B values U is undefined.  (The old b_empty
    // branch computed na and t, but then always returned HUGE_VAL via !nb.)
    for (i = 0; i < n; i++)
        if (b[i])
            break;
    if (i == n)
        return HUGE_VAL;

    // Count equal (e) and less-than (l) permutations, plus the tie
    // adjustment.
    //
    // The standard tie-corrected variance is
    //     var = na*nb/12 * ((N+1) - sum(p^3-p) / (N*(N-1)))
    // with p = a[i]+b[i] and N = na+nb.  As sum(p) == N this simplifies to
    //     var = na*nb * (N^3 - sum(p^3)) / (12*N*(N-1))
    // N^3 - sum(p^3) is accumulated incrementally: adding a bin of size p
    // to a running total S contributes (S+p)^3 - S^3 - p^3 = 3*S*p*(S+p).
    // All terms are non-negative, so there is no integer overflow, no
    // cancellation, and a single occupied bin gives exactly zero.
    int64_t e = 0, l = 0, na = 0, nb = 0;
    double d = 0;   // N^3 - sum(p^3)
    for (i = n-1; i >= 0; i--) {
        // Combinations of a[i] and b[j] for i==j
        e += (int64_t)a[i]*b[i];

        // nb is running total of b[i+1]..b[n-1].
        // Therefore a[i]*nb is the number of combinations of a[i] and b[j]
        // for all i < j.
        l += a[i]*nb; // a<b

        int64_t S = na + nb, p = (int64_t)a[i] + b[i];
        d += 3.0 * S * p * (S + p);

        na += a[i];
        nb += b[i];
    }

    if (!na)
        return HUGE_VAL;

    double N = na + nb;             // >= 2 as na, nb >= 1
    double U = l + e*0.5;           // Mann-Whitney U score
    double m = (double)na*nb / 2.0;

    // With ties adjustment
    double var2 = (double)na*nb * d / (12.0 * N * (N-1));
    // var = na*nb*(na+nb+1)/12.0; // simpler; minus tie adjustment

    if (var2 <= 0)
        return do_Z ? 0 : 1;

    if (do_Z) {
        // S.D. normalised Z-score
        //Z = (U - m - (U-m >= 0 ? 0.5 : -0.5)) / sd; // gatk method?
        return (U - m) / sqrt(var2);
    }

    // Else U score, which can be asymmetric for some data types.
    if (left_only && U > m)
        return HUGE_VAL; // one-sided, +ve bias is OK, -ve is not.

    if (na >= 8 || nb >= 8) {
        // Normal approximation, very good for na>=8 && nb>=8 and
        // reasonable if na<8 or nb<8
        return exp(-0.5*(U-m)*(U-m)/var2);
    }

    // Exact calculation; na, nb < 8 here
    if (na==1 || nb == 1)
        return mann_whitney_1947_((int)na, (int)nb, U) * sqrt(2*M_PI*var2);
    else
        return mann_whitney_1947((int)na, (int)nb, U) * sqrt(2*M_PI*var2);
}

[jkbonfield]

FYI my reproduction is:

$ cat _ref.fa
>ref
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

$ d=1300;for bcftools in ./bcftools.dev ./bcftools.PR ./bcftools;do echo === $i;perl -e 'srand(1);$d='$d';print "\@SQ\tSN:ref\tLN:1000\n";$len=10;$seq="A" x $len;$qual="I" x $len;for ($i=0;$i<$d;$i++) {$MQ=1;if (rand()>0.5) {$MQ=60;substr($seq,$len/2,1) = "C";substr($qual,$len/2, 1)="i";};print "read$i\t0\tref\t1\t$MQ\t${len}M\t*\t0\t0\t$seq\t$qual\n"}' > _.sam;  $bcftools mpileup _.sam -f _ref.fa -B -d 10000000 -a -AD 2>/dev/null|tr ';' '\012'|grep Z=;done
===
RPBZ=0
MQBZ=1.3953
BQBZ=1.74783
SCBZ=0
===
RPBZ=0
MQBZ=1.3953
BQBZ=36.0416
SCBZ=0
===
RPBZ=0
MQBZ=1.3953
BQBZ=36.0416
SCBZ=0

So I can't see a big change in MQ, but I'm unsure how this is used now without deeper digging. I can however see a significant change in base quality which occurs between depth 1292 and 1293.

At d=2097152 the BQBZ score with this PR is 1448.15. By d=2097160 it's dropped to 1.73206 while the refactored code above is unchanged. That makes sense it's scaling with (I think) approx p^2 rather than p^3.

[cindykrafft]

Thanks for digging into this. Here's the script. It takes a bcftools binary and the number of reference reads: 1300 shows the problem, 1290 doesn't.

#!/bin/sh
# Minimal reproduction: bcftools mpileup INFO/MQBZ once one MAPQ bin holds >= 1291 reads.
# Usage: sh mcve_mpileup_mqbz_overflow.sh /path/to/bcftools [NREF]   (NREF defaults to 1300; 1290 is correct)
set -e
BCF=${1:-bcftools}; NREF=${2:-1300}; NALT=40
D=$(mktemp -d); cd "$D"
# 300-bp reference "ACGTACGT..."; the site is 1-based 150 (REF C); alt reads carry T there.
awk 'BEGIN{s=""; for(i=0;i<300;i++) s=s substr("ACGT",i%4+1,1); print ">ref"; print s}' > ref.fa
printf 'ref\t300\t5\t300\t301\n' > ref.fa.fai
awk -v nref=$NREF -v nalt=$NALT 'BEGIN{
  for(i=0;i<300;i++) ref=ref substr("ACGT",i%4+1,1);
  q=""; for(i=0;i<50;i++) q=q "I";
  print "@HD\tVN:1.6\tSO:coordinate"; print "@SQ\tSN:ref\tLN:300";
  n=nref+nalt; step=int(n/nalt); a=0;
  for(i=0;i<n;i++){
    start=101+int(i*49/n); seq=substr(ref,start,50);      # sorted starts 101..149; every read covers 150
    mq=60; if(i%step==0 && a<nalt){ a++; mq=30; k=150-start+1; seq=substr(seq,1,k-1) "T" substr(seq,k+1) }
    printf "r%d\t%d\tref\t%d\t%d\t50M\t*\t0\t0\t%s\t%s\n", i, (i%2)*16, start, mq, seq, q
  }}' > reads.sam
"$BCF" mpileup -f ref.fa -B -d 100000 reads.sam 2>mpileup.err | awk '$2==150' | cut -f 2,4,5,8 | tr ';' '\n' | grep -E '^150|^DP=|MQBZ' | tr '\n' ' '; echo
# Expected: the tie-corrected Mann-Whitney U Z-score of the MAPQ bins (ref reads: bin 59, alt reads: bin 30).
awk -v na=$NREF -v nb=$NALT 'BEGIN{ N=na+nb; U=0; m=na*nb/2; T=(na^3-na)+(nb^3-nb);
  v=na*nb/12*((N+1)-T/(N*(N-1))); printf "expected MQBZ=%.4f (U=%d, mean=%d, tie-corrected var=%.2f)\n", (U-m)/sqrt(v), U, m, v }'
"$BCF" --version | head -1
cd /; rm -rf "$D"

I think your reproduction hits the same condition, just in BQ instead of MQ. In the perl one-liner $seq and $qual aren't reset inside the loop. Once the first MQ 60 read has set the C/i at position 6, every later read carries it whatever its MQ. With srand(1) only read0 and read1 are A, so at depth d one base-quality bin holds d−2 reads, and that reaches 1291 at d=1293. The MQ values split about 640/654 between MQ 1 and MQ 60, so neither MQ bin gets near 1291, and MQBZ is the same on all three builds. My numbers are the same as yours: at d=1293, BQBZ is 1.73374 on develop and 35.9444 with the PR; at d=1300 it's 1.74783 and 36.0416.

You're right that MQ isn't special. Any of the Z annotations overflows once one of its bins holds 1291 or more reads. I used MQ in the example because, with most reads at MAPQ 60, the whole reference depth tends to sit in one bin, so on real data it's the easiest one to hit.

I built your refactor on develop and checked it:

• make test: 2488 passed, 0 failed (with this PR's regression test in the tree).
• The script above: MQBZ −36.5924 at 1300 reads and −36.4555 at 1290, the same as the PR and as the tie-corrected value.
• Your reproduction: identical to the PR at d=1292, 1293 and 1300. At d=2097160 it gives BQBZ 1448.16, where the PR drops to 1.73206, as you found.
• Our harness against scipy's asymptotic tie-corrected Mann–Whitney (31 cases, including a 48-sample run with default BAQ): all agree.
• Against the PR's int64 version on about 1.9 million random histograms, some with bins of 100k–500k reads: the largest relative difference is about 1e-11, and both return HUGE_VAL in the same cases.

So I think yours is better: it removes the cube instead of moving the limit. I'm happy to replace the change in this PR with your version and keep the regression test and the NEWS entry, or to close this PR if you'd rather commit it yourself. If I update the PR, how would you like your part attributed?

---

Generated by Claude Code

[jkbonfield]

Yep I realised I had a bug in my script. I was trying to create REF/ALT SNP and it basically gave ALT only bar 1 or 2 initial reads. However correcting it gave less impactful results, so no matter the buggy one will do for validation purposes!

I think the refactored code is probably better (although harder to analyse and see what actually changed, unlike the original type change). I need to figure out if it's worth doing therefore as we don't really want to commit AI generated code without first understanding the impact of each line. A human needs to sign off on those changes basically, in order to fulfil our policies.

[jkbonfield]

So I think yours is better: it removes the cube instead of moving the limit. I'm happy to replace the change in this PR with your version and keep the regression test and the NEWS entry, or to close this PR if you'd rather commit it yourself. If I update the PR, how would you like your part attributed?

I decided in the end to just create a new PR which contains this one with the commit intact, and a follow up one. That keeps the authorship and provenance intact, simply as a way to document the progression more than anything else.

Thanks for raising the issue and providing a fix. Over to you @pd3.

[pd3]

Thank you for this. This was merged via #2598