Exact cohort distance calculations#
For four or more samples, pairwise_distances batches F/F2/D/Jaccard through
one shared clonotype vocabulary and sparse products. Up to three samples, it keeps
pairwise joins: there are at most three pairs, and a global vocabulary adds
unnecessary work for a few deep samples. Both paths aggregate each sample once.
The batch uses existing SciPy kernels, with no new native extension or fitted parameters. Pearson R retains the centered pairwise
calculation to avoid cancellation in a raw-moment covariance formula. Signed or
non-finite frequency tables also retain their existing pairwise arithmetic.
F/F2 batching uses canonical Int64 counts with effective read totals below
\(2^{53}\); other count types and larger totals retain the pairwise path.
For \(N\) samples and \(K\) fully specified clonotype keys, let \(B_{ik}\) be one when sample \(i\) contains key \(k\), and zero otherwise. Let \(W_{ik}\) be that key’s within-sample read frequency, including an explicit zero for a zero-count clonotype. Frequencies retain their original sample totals; null-key rows cannot match another sample but still contribute to those totals and the sample richness \(n_i\). Thus missing calls are not silently matched.
The implementation in overlap.cluster._exact_matrix computes
where the square root is elementwise. Entry \(C_{ij}\) counts shared keys, \(A_{ij}\) sums sample \(i\)’s frequencies on those keys, and \(H_{ij}\) sums their geometric mean frequencies. Therefore the similarities are
A zero denominator gives similarity zero. The distance transforms are \(\max(0,-\log_{10}(s+10^{-9}))\) for D/F/F2 and \(1-s\) for Jaccard, where \(s\) denotes the corresponding similarity. The zero floor prevents the logarithm’s offset from producing negative distances at similarity one, which hierarchical clustering rejects. The diagonal is set to zero and both output forms preserve sample order.
The code accumulates read counts before normalization. Let \(Q_{ik}\) be the nonnegative integer read count of fully specified key \(k\), and let \(t_i\) be the sample’s total read count, including null-key rows. A zero total uses \(t_i=1\), matching the single-pair normalization. Then
For nonnegative integer totals below \(2^{53}\), float64 shared-count accumulation is exact. This avoids accumulating millions of small normalized frequencies for F. F2 still rounds its square roots and sums. Vocabulary indices follow first occurrence, keeping accumulation order repeatable for the same inputs.
These are algebraic identities, not new estimators. Floating-point summation order can differ. For a sum of \(m\) nonnegative terms in the normal floating-point range, the usual accumulation bound is \(\gamma_m=mu/(1-mu)\) when \(mu<1\), with float64 unit roundoff \(u=2^{-53}\); products and square roots add their own rounding. A similarity perturbation from \(x\) to \(y\) changes its log distance by at most \(|x-y|/[\ln(10)(\min(x,y)+10^{-9})]\). Bitwise equality with pairwise joins is not promised. Boundary tests cover empty samples, zero counts, duplicate keys, null calls, custom keys and signed inputs.
With \(P\) aggregated sample/key rows, \(p_i\) fully specified keys in sample \(i\), and \(d_k\) samples containing key \(k\), expected fixed-length-key construction, sparse row sorting and product work is \(O(P+\sum_i p_i\log(p_i+1)+\sum_k d_k^2+N^2)\). The returned dense distances inherently require \(O(N^2)\) storage; numerical workspace is checked against the available memory allocation. Matrix and long-form outputs use native columns rather than constructing one Python float or tuple per result cell. Large all-pairs outputs still require adequate memory or a smaller sample collection.