Implement primer pair finding and evaluation logic

This change introduces the core functionality for identifying, validating, and statistically evaluating primer pairs within the corpus. It includes new data structures for primer parameters, amplicon statistics, and sequence handling, along with methods for assigning taxa to sequences. Comprehensive tests are added to ensure the accuracy of the pair evaluation against expected outputs.

OBIPRIMERS_NOTES.md

@@ -301,3 +301,59 @@ Variantes du protocole :
 
 Décision en suspens : écrire la variante 3 maintenant, ou d'abord finir le pipeline (paires et indices)
 et profiler l'ensemble pour savoir si l'appariement domine encore.
+
+## 12. Mode A, étape 4 : paires et indices (2026-11-03)
+
+Code : `taxa.go` (`Taxa`, `Corpus.AssignTaxa`), `pairs.go` (paramètres, amorces conservées, énumération
+des amplicons, contrainte du 3', comptage), `pairs_eval.go` (`FindPairs`, `PairStats`, Bc, Bs, Yule, Tm),
+test opt-in `TestPairsAgainstEcoPrimers`.
+
+Fidèle à ecoPrimers : amorce conservée = trouvée dans >= floor(|E| x 0,9) exemples ; « good » si en plus
+<= floor(|C| x 0,1) contre-exemples ; une paire a au moins une amorce good ; longueur d'amplicon
+strictement entre `-l` et `-L` ; clé de paire (amorce 1, sens, amorce 2, sens) ordonnée par mot
+canonique ; paire retenue si q > 0,9 et qfp < 0,1 (comparaisons en float32 comme le C) et Bs >= 0,6 ;
+Bc = taxons amplifiés / taxons exemples ; **Bs = taxons identifiés / taxons amplifiés par la paire**
+(et non / tous les taxons : c'est ce que fait le code, l'article est ambigu) ; un taxon est identifié si
+au moins un de ses amplicons n'est partagé avec aucun autre taxon ; tri par Yule x Bs.
+
+Choix d'implémentation (mémoire) : deux phases. Phase 1 ne garde que des compteurs par paire (séquences
+exemples/contre-exemples, longueurs) ; phase 2 ne revisite que les paires retenues (407 / 3 448) et
+n'enregistre qu'un hash 64 bits d'amplicon par occurrence au lieu de stocker les amplicons. ecoPrimers
+conserve un amplicon par occurrence de chaque paire candidate (42 682 paires pour Cyprinidae).
+
+**Contrôle contre les tables complètes d'ecoPrimers** (mêmes paramètres, `-e 3 -3 2 -l 30 -L 150 -c`) :
+
+| Jeu | Paires | Résultat |
+|---|---|---|
+| Cyprinidae | 3 448 / 3 448 | **identique sur toutes les colonnes** : amorces, good, séquences, taxons, Bc, Bs, longueurs, GC, Tm, Tm minimale |
+| Teleostei | 407 / 407 | mêmes paires ; **6 comptent 1 séquence de moins** (NC_033548) |
+
+Cause de l'écart Teleostei **établie par expérience** : l'`extractSite` d'ecoPrimers lit un `N` (code 4)
+dans un registre à décalage ; le 4 déborde sur la base voisine et `~4 & 3` vaut 3 (T), donc un `N` en 3'
+peut être accepté par accident. En émulant exactement ce registre dans le test du 3' (modification
+temporaire, annulée), les **deux tables sont reproduites à 100 %** sur toutes les colonnes. Non conservé :
+une base ambiguë ne peut pas être une correspondance 3' vérifiée. Le test fige donc l'écart (6 paires,
+-1 séquence, rien d'autre ne diffère). Une première émulation plus grossière (N lu comme A) était
+fausse : elle créait un écart dans Cyprinidae.
+
+Autres écarts connus, non testés sur ces données : fenêtres à cheval sur la jonction circulaire, que
+le C lit au-delà du tableau (comportement indéfini) alors que nous reprenons au début ; Tm minimale
+ignorant les sites contenant une base ambiguë (le C en tire des valeurs arbitraires).
+
+**Profil du pipeline complet (16 workers, M3 Max)**, le critère de l'option A :
+
+| Étape | Cyprinidae | Teleostei |
+|---|---|---|
+| lecture + sélection | 0,45 s | 0,44 s |
+| taxons | 0,00 s | 0,00 s |
+| mots stricts | 0,14 s | 1,04 s |
+| **recherche approchée** | **4,14 s** | **10,18 s** |
+| paires + indices + Tm | 1,04 s | 1,12 s |
+| total | ~5,8 s | ~12,8 s |
+
+ecoPrimers complet : 134 s et 350 s. **L'appariement approché reste 72 % à 80 % du total : la variante
+SIMD (variante 3) est justifiée par la mesure.** Mémoire : non isolée (le tas de 2,3 à 2,8 Go du test
+contient la taxonomie NCBI complète).
+
+Non fait : option `-U` (pas de multi-match), séquence de référence `-R`, sets d'amorces, sortie CSV,
+CLI. Prochaine étape : variante 3 (bit-parallèle multi-mots, SIMD).

pkg/obiprimer/pairs.go

@@ -0,0 +1,281 @@
+package obiprimer
+
+import (
+	"cmp"
+	"fmt"
+	"slices"
+	"strings"
+	"sync"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer"
+)
+
+// PairParams holds the parameters of the primer pair search, named after the
+// ecoPrimers options.
+type PairParams struct {
+	PrimerLength   int     // -O, word size
+	MinLength      int     // -l: the amplicon (primers excluded) must be longer than this
+	MaxLength      int     // -L: and shorter than this
+	ThreePrime     int     // -3: number of 3' bases that must match exactly (0 = no constraint)
+	Sensitivity    float64 // -s: a primer must hit more than this fraction of the examples (0.9)
+	FalsePositive  float64 // -x: and less than this fraction of the counterexamples (0.1)
+	MinSpecificity float64 // -T: pairs with a lower Bs are dropped (0.6)
+	Circular       bool    // -c
+}
+
+func (p PairParams) validate() error {
+	switch {
+	case p.PrimerLength < 2 || p.PrimerLength > 31:
+		return fmt.Errorf("obiprimer: primer length must be between 2 and 31, got %d", p.PrimerLength)
+	case p.MaxLength <= p.MinLength+1:
+		return fmt.Errorf("obiprimer: empty amplicon length range ]%d,%d[", p.MinLength, p.MaxLength)
+	case p.ThreePrime < 0 || p.ThreePrime > p.PrimerLength:
+		return fmt.Errorf("obiprimer: 3' strict bases must be in [0,%d], got %d", p.PrimerLength, p.ThreePrime)
+	case p.Sensitivity <= 0 || p.Sensitivity > 1 || p.FalsePositive < 0 || p.FalsePositive > 1:
+		return fmt.Errorf("obiprimer: quorums must be in ]0,1] and [0,1]")
+	}
+	return nil
+}
+
+// Primer is a conserved word: one found, with at most the allowed errors, in
+// enough example sequences to be a candidate primer.
+type Primer struct {
+	Word            uint64 // canonical word
+	Examples        int    // example sequences where it is found (either strand)
+	CounterExamples int
+	// Good means conserved in the examples and rare enough in the
+	// counterexamples; a pair needs at least one good primer.
+	Good bool
+}
+
+// seqMatch is one occurrence of a primer in one sequence.
+type seqMatch struct {
+	primer  int32
+	pos     uint32
+	reverse bool
+}
+
+// primerSet is the conserved words and, for each corpus sequence, their
+// occurrences sorted by position.
+type primerSet struct {
+	primers []Primer
+	perSeq  [][]seqMatch
+}
+
+// quorumCount is floor(n*fraction), robust to the binary form of the fraction.
+func quorumCount(n int, fraction float64) int { return MinSequences(n, fraction) }
+
+// selectPrimers keeps the words found in at least floor(examples*Sensitivity)
+// example sequences and marks as good those found in at most
+// floor(counterexamples*FalsePositive) counterexamples.
+func selectPrimers(c *Corpus, wm []WordMatches, p PairParams) *primerSet {
+	counts := c.Counts()
+	inQ := quorumCount(counts.Examples, p.Sensitivity)
+	outQ := quorumCount(counts.CounterExamples, p.FalsePositive)
+
+	ps := &primerSet{perSeq: make([][]seqMatch, len(c.Seqs))}
+	for _, w := range wm {
+		ex, co := w.SequencesMatched(c)
+		if ex < inQ {
+			continue
+		}
+		id := int32(len(ps.primers))
+		ps.primers = append(ps.primers, Primer{w.Word, ex, co, co <= outQ})
+		for _, m := range w.Matches {
+			ps.perSeq[m.Seq] = append(ps.perSeq[m.Seq], seqMatch{id, m.Pos, m.Reverse})
+		}
+	}
+	for _, ms := range ps.perSeq {
+		slices.SortStableFunc(ms, func(a, b seqMatch) int { return cmp.Compare(a.pos, b.pos) })
+	}
+	return ps
+}
+
+// pairKey identifies a pair: the two primers, ordered by canonical word, and
+// for each whether the primer is the canonical word as is (true) or its reverse
+// complement (false). It is the key of ecoPrimers' pair tree.
+type pairKey struct {
+	p1, p2 int32
+	a1, a2 bool
+}
+
+// occurrence is one amplicon: a pair of matches in a sequence.
+type occurrence struct {
+	key      pairKey
+	swapped  bool // true when the upstream match is the second primer of the key
+	up, down seqMatch
+	length   int // amplicon length, primers excluded
+}
+
+// forEachOccurrence calls fn for every amplicon of sequence s, applying the
+// length, good-primer and 3' constraints of ecoPrimers.
+func (ps *primerSet) forEachOccurrence(c *Corpus, s int, p PairParams, fn func(occurrence)) {
+	ms := ps.perSeq[s]
+	lp := p.PrimerLength
+	var raw []byte
+	if p.ThreePrime > 0 {
+		raw = c.Seqs[s].Sequence()
+	}
+
+	for i := range ms {
+		for j := i + 1; j < len(ms); j++ {
+			dpos := int(ms[j].pos) - int(ms[i].pos)
+			if dpos <= lp {
+				continue
+			}
+			length := dpos - lp
+			if length >= p.MaxLength {
+				break
+			}
+			up, down := ms[i], ms[j]
+			if length <= p.MinLength || !(ps.primers[up.primer].Good || ps.primers[down.primer].Good) {
+				continue
+			}
+
+			// The upstream match is read as is (primer 1 of the pair when its word
+			// is smaller), the downstream one as the reverse primer.
+			key := pairKey{up.primer, down.primer, !up.reverse, down.reverse}
+			swapped := !(ps.primers[key.p2].Word > ps.primers[key.p1].Word)
+			if swapped {
+				key = pairKey{key.p2, key.p1, key.a2, key.a1}
+			}
+
+			if p.ThreePrime > 0 && !ps.threePrimeMatches(raw, up, down, p) {
+				continue
+			}
+			fn(occurrence{key, swapped, up, down, length})
+		}
+	}
+}
+
+// wordBase returns base t (0 = 5') of a 2-bit packed word of size lp.
+func wordBase(word uint64, lp, t int) byte { return byte(word>>(2*uint(lp-1-t))) & 3 }
+
+// threePrimeMatches reports whether the last p.ThreePrime bases of the forward
+// primer (as found at the upstream match) and of the reverse primer (reverse
+// complement of the downstream window) equal the sequence, i.e. have no
+// mismatch at their 3' end.
+func (ps *primerSet) threePrimeMatches(raw []byte, up, down seqMatch, p PairParams) bool {
+	lp, n := p.PrimerLength, len(raw)
+	at := func(pos int) byte {
+		if p.Circular {
+			pos %= n
+		}
+		return baseCode[raw[pos]]
+	}
+
+	fw := ps.primers[up.primer].Word
+	rv := ps.primers[down.primer].Word
+	for t := lp - p.ThreePrime; t < lp; t++ {
+		// forward primer: the word when found direct, its reverse complement otherwise
+		var f byte
+		if !up.reverse {
+			f = wordBase(fw, lp, t)
+		} else {
+			f = 3 - wordBase(fw, lp, lp-1-t)
+		}
+		if s := at(int(up.pos) + t); s > 3 || s != f {
+			return false
+		}
+
+		// reverse primer: the word when found reverse, its reverse complement otherwise
+		var r byte
+		if down.reverse {
+			r = wordBase(rv, lp, t)
+		} else {
+			r = 3 - wordBase(rv, lp, lp-1-t)
+		}
+		if s := at(int(down.pos) + lp - 1 - t); s > 3 || 3-s != r {
+			return false
+		}
+	}
+	return true
+}
+
+// pairAcc accumulates, over the sequences, the counts ecoPrimers keeps while it
+// builds the pairs.
+type pairAcc struct {
+	examples, counters uint32
+	last               int32 // 1 + last sequence counted, to count each sequence once
+	amplicons          uint32
+	sumLen             uint64
+	minLen, maxLen     uint32
+}
+
+// countPairs enumerates the amplicons of every sequence, in parallel, and
+// returns for each pair the sequences it amplifies and the length statistics
+// of its example amplicons.
+func (ps *primerSet) countPairs(c *Corpus, p PairParams) map[pairKey]*pairAcc {
+	workers := obidefault.ParallelWorkers()
+	partial := make([]map[pairKey]*pairAcc, workers)
+
+	jobs := make(chan int)
+	var wg sync.WaitGroup
+	for w := range workers {
+		wg.Add(1)
+		partial[w] = map[pairKey]*pairAcc{}
+		go func() {
+			defer wg.Done()
+			acc := partial[w]
+			for s := range jobs {
+				isExample := c.Roles[s] == Example
+				ps.forEachOccurrence(c, s, p, func(o occurrence) {
+					a := acc[o.key]
+					if a == nil {
+						a = &pairAcc{minLen: ^uint32(0)}
+						acc[o.key] = a
+					}
+					if a.last != int32(s)+1 {
+						a.last = int32(s) + 1
+						if isExample {
+							a.examples++
+						} else {
+							a.counters++
+						}
+					}
+					if isExample {
+						l := uint32(o.length)
+						a.amplicons++
+						a.sumLen += uint64(l)
+						a.minLen = min(a.minLen, l)
+						a.maxLen = max(a.maxLen, l)
+					}
+				})
+			}
+		}()
+	}
+	for s := range c.Seqs {
+		jobs <- s
+	}
+	close(jobs)
+	wg.Wait()
+
+	// A sequence is handled by one worker only, so the sequence counts add up.
+	total := partial[0]
+	for _, m := range partial[1:] {
+		for k, a := range m {
+			t := total[k]
+			if t == nil {
+				total[k] = a
+				continue
+			}
+			t.examples += a.examples
+			t.counters += a.counters
+			t.amplicons += a.amplicons
+			t.sumLen += a.sumLen
+			t.minLen = min(t.minLen, a.minLen)
+			t.maxLen = max(t.maxLen, a.maxLen)
+		}
+	}
+	return total
+}
+
+// orientedWord returns the upper-case primer string 5'->3': the canonical word
+// when direct, its reverse complement otherwise.
+func orientedWord(word uint64, lp int, direct bool) string {
+	if !direct {
+		word = obikmer.ReverseComplement(word, lp)
+	}
+	return strings.ToUpper(string(obikmer.DecodeKmer(word, lp, nil)))
+}

pkg/obiprimer/pairs_eval.go

@@ -0,0 +1,360 @@
+package obiprimer
+
+import (
+	"cmp"
+	"fmt"
+	"math"
+	"slices"
+	"sync"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obithermo"
+)
+
+// PairStats describes one primer pair and its quality indices. Primers are
+// written 5'->3' and ordered so that Primer1 < Primer2 (as in ecoPrimers).
+type PairStats struct {
+	Primer1, Primer2 string
+	Good1, Good2     bool // primer conserved in the examples and rare in the counterexamples
+
+	Examples        int // example sequences amplified
+	CounterExamples int
+	Yule            float32 // Examples/|E| - CounterExamples/|C|
+
+	ExampleTaxa    int     // taxa amplified among the examples
+	CounterTaxa    int     // taxa amplified among the counterexamples
+	Bc             float32 // coverage: ExampleTaxa over all the example taxa
+	Identified     int     // example taxa with an amplicon no other taxon shares
+	Bs             float32 // specificity: Identified over ExampleTaxa
+	MinLen, MaxLen int
+	MeanLen        float64
+	Tm1, Tm2       float64 // melting temperature of the primers against their perfect target
+	MinTm1, MinTm2 float64 // lowest melting temperature over the example amplicons
+	GC1, GC2       int
+}
+
+// Score is the ranking criterion of ecoPrimers.
+func (p PairStats) Score() float32 { return p.Yule * p.Bs }
+
+// ampRecord is an amplicon of a retained pair, in one sequence.
+type ampRecord struct {
+	pair     int32
+	seq      int32
+	hash     uint64
+	tm1, tm2 float32 // NaN when unavailable
+}
+
+// ampHash hashes the amplicon between the two primers read from the side of
+// the first primer of the pair (reverse complemented when the upstream match is
+// the second primer), so that identical amplicons give the same value
+// whatever their strand. Non ACGT bases hash as a fifth symbol.
+func ampHash(raw []byte, from, to int, reverse bool) uint64 {
+	const prime = 1099511628211
+	h := uint64(14695981039346656037) ^ uint64(to-from)
+	h *= prime
+	if !reverse {
+		for _, b := range raw[from:to] {
+			h = (h ^ uint64(baseCode[b])) * prime
+		}
+		return h
+	}
+	for i := to - 1; i >= from; i-- {
+		c := baseCode[raw[i]]
+		if c < 4 {
+			c = 3 - c
+		}
+		h = (h ^ uint64(c)) * prime
+	}
+	return h
+}
+
+// FindPairs searches the primer pairs from the approximate occurrences of the
+// strict words (see MatchApat) and evaluates them.
+//
+// A pair is a pair of conserved primers, at least one of them good, that
+// delimit an amplicon of the allowed length. It is retained when it amplifies
+// more than Sensitivity of the examples and less than FalsePositive of the
+// counterexamples, and its specificity Bs reaches MinSpecificity. The result
+// is sorted by decreasing Yule*Bs.
+//
+// model computes the melting temperatures; nil leaves them NaN.
+func FindPairs(c *Corpus, taxa *Taxa, wm []WordMatches, p PairParams, model *obithermo.Model) ([]PairStats, error) {
+	if err := p.validate(); err != nil {
+		return nil, err
+	}
+	if taxa == nil || len(taxa.Of) != len(c.Seqs) {
+		return nil, fmt.Errorf("obiprimer: taxa do not match the corpus")
+	}
+	ps := selectPrimers(c, wm, p)
+	if len(ps.primers) == 0 {
+		return nil, nil
+	}
+
+	// Phase 1: which pairs amplify enough examples and few enough counterexamples.
+	counts := c.Counts()
+	acc := ps.countPairs(c, p)
+	sens, fpq := float32(p.Sensitivity), float32(p.FalsePositive)
+	var kept []pairKey
+	for k, a := range acc {
+		q := float32(a.examples) / float32(counts.Examples)
+		qfp := float32(0)
+		if counts.CounterExamples > 0 {
+			qfp = float32(a.counters) / float32(counts.CounterExamples)
+		}
+		if q > sens && qfp < fpq {
+			kept = append(kept, k)
+		}
+	}
+	slices.SortFunc(kept, func(a, b pairKey) int {
+		return cmp.Or(cmp.Compare(a.p1, b.p1), cmp.Compare(a.p2, b.p2),
+			boolCmp(a.a1, b.a1), boolCmp(a.a2, b.a2))
+	})
+	index := make(map[pairKey]int32, len(kept))
+	for i, k := range kept {
+		index[k] = int32(i)
+	}
+
+	// Phase 2: the amplicons of the retained pairs only.
+	records := ps.collectAmplicons(c, p, index, model)
+	byPair := make([][]ampRecord, len(kept))
+	for _, r := range records {
+		byPair[r.pair] = append(byPair[r.pair], r)
+	}
+
+	var out []PairStats
+	for i, k := range kept {
+		st, ok := ps.evaluate(c, taxa, p, k, acc[k], byPair[i], model)
+		if ok {
+			out = append(out, st)
+		}
+	}
+	slices.SortFunc(out, func(a, b PairStats) int {
+		return cmp.Or(-cmp.Compare(a.Score(), b.Score()), cmp.Compare(a.Primer1, b.Primer1), cmp.Compare(a.Primer2, b.Primer2))
+	})
+	return out, nil
+}
+
+func boolCmp(a, b bool) int {
+	switch {
+	case a == b:
+		return 0
+	case !a:
+		return -1
+	}
+	return 1
+}
+
+// collectAmplicons gathers, in parallel, one record per amplicon of a retained pair.
+func (ps *primerSet) collectAmplicons(c *Corpus, p PairParams, index map[pairKey]int32, model *obithermo.Model) []ampRecord {
+	lp := p.PrimerLength
+	workers := obidefault.ParallelWorkers()
+	parts := make([][]ampRecord, workers)
+
+	// oriented[primer][0] is the reverse complement of the word, [1] the word itself.
+	oriented := make([][2][]byte, len(ps.primers))
+	for i, pr := range ps.primers {
+		oriented[i][0] = []byte(orientedWord(pr.Word, lp, false))
+		oriented[i][1] = []byte(orientedWord(pr.Word, lp, true))
+	}
+	nan := float32(math.NaN())
+
+	jobs := make(chan int)
+	var wg sync.WaitGroup
+	for w := range workers {
+		wg.Add(1)
+		go func() {
+			defer wg.Done()
+			var site1, site2 [32]byte
+			for s := range jobs {
+				raw := c.Seqs[s].Sequence()
+				n := len(raw)
+				isExample := c.Roles[s] == Example
+				ps.forEachOccurrence(c, s, p, func(o occurrence) {
+					pi, ok := index[o.key]
+					if !ok {
+						return
+					}
+					rec := ampRecord{pair: pi, seq: int32(s), tm1: nan, tm2: nan}
+					if isExample {
+						rec.hash = ampHash(raw, int(o.up.pos)+lp, int(o.down.pos), o.swapped)
+						if model != nil {
+							fp := oriented[o.up.primer][b2i(!o.up.reverse)]
+							rp := oriented[o.down.primer][b2i(o.down.reverse)]
+							fsite, rsite := site1[:lp], site2[:lp]
+							okF := windowAsIs(fsite, raw, int(o.up.pos), n, p.Circular)
+							okR := windowRevComp(rsite, raw, int(o.down.pos), n, p.Circular)
+							tmF, tmR := nan, nan
+							if okF {
+								if t, err := model.Tm(fp, fsite); err == nil {
+									tmF = float32(t)
+								}
+							}
+							if okR {
+								if t, err := model.Tm(rp, rsite); err == nil {
+									tmR = float32(t)
+								}
+							}
+							if o.swapped {
+								rec.tm1, rec.tm2 = tmR, tmF
+							} else {
+								rec.tm1, rec.tm2 = tmF, tmR
+							}
+						}
+					}
+					parts[w] = append(parts[w], rec)
+				})
+			}
+		}()
+	}
+	for s := range c.Seqs {
+		jobs <- s
+	}
+	close(jobs)
+	wg.Wait()
+
+	var all []ampRecord
+	for _, part := range parts {
+		all = append(all, part...)
+	}
+	return all
+}
+
+func b2i(b bool) int {
+	if b {
+		return 1
+	}
+	return 0
+}
+
+const upper = "ACGT"
+
+// windowAsIs copies the len(dst) bases at pos as upper-case ACGT; it reports
+// false if one of them is not A, C, G or T.
+func windowAsIs(dst, raw []byte, pos, n int, circular bool) bool {
+	for t := range dst {
+		i := pos + t
+		if circular {
+			i %= n
+		}
+		c := baseCode[raw[i]]
+		if c > 3 {
+			return false
+		}
+		dst[t] = upper[c]
+	}
+	return true
+}
+
+// windowRevComp writes the reverse complement of the len(dst) bases at pos.
+func windowRevComp(dst, raw []byte, pos, n int, circular bool) bool {
+	l := len(dst)
+	for t := range dst {
+		i := pos + l - 1 - t
+		if circular {
+			i %= n
+		}
+		c := baseCode[raw[i]]
+		if c > 3 {
+			return false
+		}
+		dst[t] = upper[3-c]
+	}
+	return true
+}
+
+// evaluate computes the indices of one retained pair. It reports false when the
+// pair must be dropped (no taxon amplified, or specificity under the threshold).
+func (ps *primerSet) evaluate(c *Corpus, taxa *Taxa, p PairParams, k pairKey, a *pairAcc,
+	recs []ampRecord, model *obithermo.Model) (PairStats, bool) {
+
+	exTaxa := map[int32]struct{}{}
+	coTaxa := map[int32]struct{}{}
+	ampTaxon := map[uint64]int32{} // amplicon -> its taxon, or -2 when shared
+	for _, r := range recs {
+		t := taxa.Of[r.seq]
+		if t < 0 {
+			continue
+		}
+		if c.Roles[r.seq] != Example {
+			coTaxa[t] = struct{}{}
+			continue
+		}
+		exTaxa[t] = struct{}{}
+		if prev, seen := ampTaxon[r.hash]; !seen {
+			ampTaxon[r.hash] = t
+		} else if prev != t {
+			ampTaxon[r.hash] = -2
+		}
+	}
+	if len(exTaxa) == 0 {
+		return PairStats{}, false
+	}
+	identified := map[int32]struct{}{}
+	for _, t := range ampTaxon {
+		if t >= 0 {
+			identified[t] = struct{}{}
+		}
+	}
+
+	counts := c.Counts()
+	st := PairStats{
+		Examples:        int(a.examples),
+		CounterExamples: int(a.counters),
+		ExampleTaxa:     len(exTaxa),
+		CounterTaxa:     len(coTaxa),
+		Identified:      len(identified),
+		MinLen:          int(a.minLen),
+		MaxLen:          int(a.maxLen),
+		MeanLen:         float64(a.sumLen) / float64(a.amplicons),
+	}
+	q := float32(a.examples) / float32(counts.Examples)
+	qfp := float32(0)
+	if counts.CounterExamples > 0 {
+		qfp = float32(a.counters) / float32(counts.CounterExamples)
+	}
+	st.Yule = q - qfp
+	st.Bc = float32(st.ExampleTaxa) / float32(taxa.ExampleTaxa)
+	st.Bs = float32(st.Identified) / float32(st.ExampleTaxa)
+	if st.Bs < float32(p.MinSpecificity) {
+		return PairStats{}, false
+	}
+
+	lp := p.PrimerLength
+	st.Primer1 = orientedWord(ps.primers[k.p1].Word, lp, k.a1)
+	st.Primer2 = orientedWord(ps.primers[k.p2].Word, lp, k.a2)
+	st.Good1, st.Good2 = ps.primers[k.p1].Good, ps.primers[k.p2].Good
+	st.GC1, st.GC2 = obithermo.GCCount([]byte(st.Primer1)), obithermo.GCCount([]byte(st.Primer2))
+
+	st.Tm1, st.Tm2, st.MinTm1, st.MinTm2 = math.NaN(), math.NaN(), math.NaN(), math.NaN()
+	if model != nil {
+		st.Tm1, _ = model.SelfTm([]byte(st.Primer1))
+		st.Tm2, _ = model.SelfTm([]byte(st.Primer2))
+		m1, m2 := float32(math.Inf(1)), float32(math.Inf(1))
+		for _, r := range recs {
+			if c.Roles[r.seq] != Example {
+				continue
+			}
+			if r.tm1 == r.tm1 { // not NaN
+				m1 = min(m1, r.tm1)
+			}
+			if r.tm2 == r.tm2 {
+				m2 = min(m2, r.tm2)
+			}
+		}
+		if !math.IsInf(float64(m1), 1) {
+			st.MinTm1 = float64(m1)
+		}
+		if !math.IsInf(float64(m2), 1) {
+			st.MinTm2 = float64(m2)
+		}
+	}
+
+	if st.Primer2 < st.Primer1 {
+		st.Primer1, st.Primer2 = st.Primer2, st.Primer1
+		st.Good1, st.Good2 = st.Good2, st.Good1
+		st.GC1, st.GC2 = st.GC2, st.GC1
+		st.Tm1, st.Tm2 = st.Tm2, st.Tm1
+		st.MinTm1, st.MinTm2 = st.MinTm2, st.MinTm1
+	}
+	return st, true
+}

pkg/obiprimer/pairs_test.go

@@ -0,0 +1,237 @@
+package obiprimer
+
+import (
+	"bufio"
+	"fmt"
+	"math"
+	"os"
+	"path/filepath"
+	"runtime"
+	"strconv"
+	"strings"
+	"testing"
+	"time"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiformats"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obithermo"
+)
+
+// ecoRow is one result line of ecoPrimers.
+type ecoRow struct {
+	p1, p2         string
+	tm1, mtm1      float64
+	tm2, mtm2      float64
+	gc1, gc2       int
+	good           string
+	in, out        int
+	yule           float64
+	intaxa, outtax int
+	bc             float64
+	well           int
+	bs             float64
+	mind, maxd     int
+	avg            float64
+}
+
+func readEcoRows(t *testing.T, path string) []ecoRow {
+	t.Helper()
+	f, err := os.Open(path)
+	if err != nil {
+		t.Fatal(err)
+	}
+	defer f.Close()
+	var rows []ecoRow
+	sc := bufio.NewScanner(f)
+	sc.Buffer(make([]byte, 1<<20), 1<<20)
+	for sc.Scan() {
+		line := sc.Text()
+		if strings.HasPrefix(line, "#") || strings.TrimSpace(line) == "" {
+			continue
+		}
+		f := strings.Split(line, "\t")
+		fl := func(i int) float64 {
+			v, err := strconv.ParseFloat(strings.TrimSpace(f[i]), 64)
+			if err != nil {
+				t.Fatalf("field %d of %q: %v", i, line, err)
+			}
+			return v
+		}
+		in := func(i int) int { return int(fl(i)) }
+		rows = append(rows, ecoRow{
+			p1: f[1], p2: f[2], tm1: fl(3), mtm1: fl(4), tm2: fl(5), mtm2: fl(6),
+			gc1: in(7), gc2: in(8), good: f[9], in: in(10), out: in(11), yule: fl(12),
+			intaxa: in(13), outtax: in(14), bc: fl(15), well: in(16), bs: fl(17),
+			mind: in(18), maxd: in(19), avg: fl(20),
+		})
+	}
+	if err := sc.Err(); err != nil {
+		t.Fatal(err)
+	}
+	return rows
+}
+
+// Opt-in end-to-end comparison with the complete ecoPrimers result tables of
+// the tutorial: same parameters (-e 3 -3 2 -l 30 -L 150 -c, -r taxid).
+//
+// Cyprinidae is reproduced exactly: its 3448 pairs agree on every column
+// (primers, good flags, sequence and taxon counts, Bc, Bs, lengths, GC, Tm and
+// lowest Tm). Teleostei finds the same 407 pairs, and 6 of them count one
+// sequence less than ecoPrimers (NC_033548, whose downstream site is
+// NGGGGAGGCAAGTCGTAA). The 3' check of ecoPrimers reads an N through the
+// shift register of extractSite, where the code 4 of N spills on the next base
+// and ~4&3 reads as T, so an N can be accepted by accident. Emulating that
+// register exactly in the 3' check was verified to reproduce both tables on
+// every column; it is not kept: an ambiguous base cannot be a verified 3' match.
+func TestPairsAgainstEcoPrimers(t *testing.T) {
+	dir := os.Getenv("OBIPRIMERS_TUTORIAL")
+	if dir == "" {
+		t.Skip("OBIPRIMERS_TUTORIAL not set")
+	}
+	taxonomy, err := obiformats.LoadTaxonomy(filepath.Join(dir, "ncbitaxo_20250211.tgz"), false, false)
+	if err != nil {
+		t.Fatal(err)
+	}
+
+	for _, tc := range []struct {
+		name, taxid string
+		species     int
+		divergent   int
+	}{{"Cyprinidae", "7953", 316, 0}, {"Teleostei", "32443", 3876, 6}} {
+		t.Run(tc.name, func(t *testing.T) {
+			var stages []string
+			tick := time.Now()
+			lap := func(name string) {
+				stages = append(stages, fmt.Sprintf("%s %s", name, time.Since(tick).Round(time.Millisecond)))
+				tick = time.Now()
+			}
+			restrict, err := BelongingTo(taxonomy, []string{tc.taxid})
+			if err != nil {
+				t.Fatal(err)
+			}
+			seqs, err := obiformats.ReadSequencesFromFile(filepath.Join(dir, "mito.vert.fasta.gz"))
+			if err != nil {
+				t.Fatal(err)
+			}
+			corpus, err := LoadCorpus(seqs, Selection{Restrict: restrict})
+			if err != nil {
+				t.Fatal(err)
+			}
+			lap("read+select")
+			taxa, err := corpus.AssignTaxa(taxonomy, "species")
+			if err != nil {
+				t.Fatal(err)
+			}
+			lap("taxa")
+			if taxa.ExampleTaxa != tc.species {
+				t.Errorf("%d example species, ecoPrimers reports %d", taxa.ExampleTaxa, tc.species)
+			}
+
+			sw, err := FindStrictWords(corpus, WordSearch{PrimerLength: 18, Quorum: 0.7, Circular: true, TmpDir: t.TempDir()})
+			if err != nil {
+				t.Fatal(err)
+			}
+			lap("strict words")
+			words := make([]uint64, len(sw.Words))
+			for i, w := range sw.Words {
+				words[i] = w.Kmer
+			}
+			wm, err := matchApat(corpus, words, MatchParams{K: 18, MaxErrors: 3, Circular: true})
+			if err != nil {
+				t.Fatal(err)
+			}
+			lap("approximate search")
+			got, err := FindPairs(corpus, taxa, wm, PairParams{
+				PrimerLength: 18, MinLength: 30, MaxLength: 150, ThreePrime: 2,
+				Sensitivity: 0.9, FalsePositive: 0.1, MinSpecificity: 0.6, Circular: true,
+			}, obithermo.DefaultModel())
+			if err != nil {
+				t.Fatal(err)
+			}
+
+			lap("pairs+indices+Tm")
+			var ms runtime.MemStats
+			runtime.ReadMemStats(&ms)
+			t.Logf("stages: %s; Go heap in use at the end %.0f MB (includes the taxonomy)", strings.Join(stages, ", "), float64(ms.HeapAlloc)/1e6)
+
+			want := readEcoRows(t, filepath.Join(dir, tc.name+".ecoprimers"))
+			t.Logf("%d pairs, ecoPrimers %d", len(got), len(want))
+
+			byKey := map[string]PairStats{}
+			for _, g := range got {
+				byKey[g.Primer1+"/"+g.Primer2] = g
+			}
+			bad := map[string]int{}
+			missing, divergent := 0, 0
+			note := func(field string, w ecoRow, format string, args ...any) {
+				if bad[field]++; bad[field] <= 3 {
+					t.Errorf("%s/%s %s: "+format, append([]any{w.p1, w.p2, field}, args...)...)
+				}
+			}
+			near := func(a, b, tol float64) bool { return math.Abs(a-b) <= tol }
+			for _, w := range want {
+				g, ok := byKey[w.p1+"/"+w.p2]
+				if !ok {
+					if missing++; missing <= 5 {
+						t.Errorf("ecoPrimers pair %s/%s not found", w.p1, w.p2)
+					}
+					continue
+				}
+				delete(byKey, w.p1+"/"+w.p2)
+				good := string("bG"[b2i(g.Good1)]) + string("bG"[b2i(g.Good2)])
+				if good != w.good {
+					note("good", w, "got %s want %s", good, w.good)
+				}
+				// Known divergence, see the comment of the test: ecoPrimers accepts one
+				// amplicon (NC_033548) whose 3' primer site ends with an N. We reject it,
+				// so the counts can be exactly one lower and nothing else may differ.
+				lower := func(got, want int) bool { return got == want || got == want-1 }
+				diverges := g.Examples != w.in || g.ExampleTaxa != w.intaxa || g.Identified != w.well
+				if diverges {
+					divergent++
+				}
+				if !lower(g.Examples, w.in) || g.CounterExamples != w.out {
+					note("sequences", w, "got %d/%d want %d/%d", g.Examples, g.CounterExamples, w.in, w.out)
+				}
+				if !lower(g.ExampleTaxa, w.intaxa) || g.CounterTaxa != w.outtax {
+					note("taxa", w, "got %d/%d want %d/%d", g.ExampleTaxa, g.CounterTaxa, w.intaxa, w.outtax)
+				}
+				if !lower(g.Identified, w.well) {
+					note("identified", w, "got %d want %d", g.Identified, w.well)
+				}
+				tol := 0.0006
+				if diverges {
+					tol = 0.0021
+				}
+				if !near(float64(g.Yule), w.yule, tol) || !near(float64(g.Bc), w.bc, tol) || !near(float64(g.Bs), w.bs, tol) {
+					note("indices", w, "yule/bc/bs got %.4f/%.4f/%.4f want %.3f/%.3f/%.3f", g.Yule, g.Bc, g.Bs, w.yule, w.bc, w.bs)
+				}
+				if g.MinLen != w.mind || g.MaxLen != w.maxd || !near(g.MeanLen, w.avg, 0.006) {
+					note("lengths", w, "got %d/%d/%.2f want %d/%d/%.2f", g.MinLen, g.MaxLen, g.MeanLen, w.mind, w.maxd, w.avg)
+				}
+				if g.GC1 != w.gc1 || g.GC2 != w.gc2 {
+					note("gc", w, "got %d/%d want %d/%d", g.GC1, g.GC2, w.gc1, w.gc2)
+				}
+				// ecoPrimers prints Kelvin-273.0; obithermo uses 273.15
+				if !near(g.Tm1+0.15, w.tm1, 0.051) || !near(g.Tm2+0.15, w.tm2, 0.051) {
+					note("tm", w, "got %.2f/%.2f want %.1f/%.1f", g.Tm1+0.15, g.Tm2+0.15, w.tm1, w.tm2)
+				}
+				if !near(g.MinTm1+0.15, w.mtm1, 0.051) || !near(g.MinTm2+0.15, w.mtm2, 0.051) {
+					note("mintm", w, "got %.2f/%.2f want %.1f/%.1f", g.MinTm1+0.15, g.MinTm2+0.15, w.mtm1, w.mtm2)
+				}
+			}
+			t.Logf("missing from ours: %d; extra in ours: %d; mismatching fields: %v; pairs one sequence below ecoPrimers: %d",
+				missing, len(byKey), bad, divergent)
+			if divergent != tc.divergent {
+				t.Errorf("%d pairs differ from ecoPrimers by one sequence, expected %d", divergent, tc.divergent)
+			}
+			for k := range byKey {
+				if len(byKey) > 0 {
+					t.Errorf("extra pair %s", k)
+				}
+				break
+			}
+		})
+	}
+	_ = obikmer.DecodeKmer
+}

pkg/obiprimer/taxa.go

@@ -0,0 +1,69 @@
+package obiprimer
+
+import (
+	"fmt"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitax"
+)
+
+// Taxa gives every corpus sequence its taxon at one rank (the species by
+// default in ecoPrimers), as a dense index, and the number of distinct taxa
+// among the examples and among the counterexamples.
+//
+// A sequence whose taxon has no ancestor at that rank (or no taxon at all) has
+// the index -1; it is never counted as a taxon.
+type Taxa struct {
+	Rank            string
+	Of              []int32 // Of[i] is the taxon index of Corpus.Seqs[i], or -1
+	ExampleTaxa     int
+	CounterExamples int // distinct taxa among the counterexamples
+	nodes           []*obitax.Taxon
+}
+
+// Name returns the taxon with the given index, for reporting.
+func (t *Taxa) Name(i int32) *obitax.Taxon { return t.nodes[i] }
+
+// AssignTaxa determines the taxon of every sequence of the corpus at rank.
+func (c *Corpus) AssignTaxa(taxonomy *obitax.Taxonomy, rank string) (*Taxa, error) {
+	if taxonomy == nil {
+		return nil, fmt.Errorf("obiprimer: a taxonomy is required to compute taxon statistics")
+	}
+	ranks := taxonomy.RankList()
+	known := false
+	for _, r := range ranks {
+		if r == rank {
+			known = true
+			break
+		}
+	}
+	if !known {
+		return nil, fmt.Errorf("obiprimer: unknown taxonomic rank %q", rank)
+	}
+
+	t := &Taxa{Rank: rank, Of: make([]int32, len(c.Seqs))}
+	index := map[*obitax.TaxNode]int32{}
+	exSeen := map[int32]struct{}{}
+	coSeen := map[int32]struct{}{}
+
+	for i, seq := range c.Seqs {
+		t.Of[i] = -1
+		taxon := seq.Taxon(taxonomy).TaxonAtRank(rank)
+		if taxon == nil {
+			continue
+		}
+		id, ok := index[taxon.Node]
+		if !ok {
+			id = int32(len(t.nodes))
+			index[taxon.Node] = id
+			t.nodes = append(t.nodes, taxon)
+		}
+		t.Of[i] = id
+		if c.Roles[i] == Example {
+			exSeen[id] = struct{}{}
+		} else {
+			coSeen[id] = struct{}{}
+		}
+	}
+	t.ExampleTaxa, t.CounterExamples = len(exSeen), len(coSeen)
+	return t, nil
+}