From 56f5fe33d726b493a6e0bcfdbe095068421db5f1 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Tue, 3 Nov 2026 14:18:39 +1100 Subject: [PATCH] Implement strict word finding and sequence filtering This introduces new structures for defining sequence roles and constraints, core logic for strict word identification, and comprehensive unit tests to ensure accuracy and correct behavior. --- OBIPRIMERS_NOTES.md | 38 ++++++ pkg/obiprimer/roles.go | 76 ++++++++++++ pkg/obiprimer/words.go | 118 +++++++++++++++++++ pkg/obiprimer/words_test.go | 223 ++++++++++++++++++++++++++++++++++++ 4 files changed, 455 insertions(+) create mode 100644 pkg/obiprimer/roles.go create mode 100644 pkg/obiprimer/words.go create mode 100644 pkg/obiprimer/words_test.go diff --git a/OBIPRIMERS_NOTES.md b/OBIPRIMERS_NOTES.md index a1d88a0..450a692 100644 --- a/OBIPRIMERS_NOTES.md +++ b/OBIPRIMERS_NOTES.md @@ -211,3 +211,41 @@ Un bug trouvé par ces tests dans ma boucle de réduction (compteur de séquence Contrôle de bout en bout prévu (étape suivante) : `Strict primer count` d'ecoPrimers = **1 208** (Cyprinidae) et **255** (Teleostei) avec `-e 3 -3 2 -l 30 -L 150 -c`, k = 18, quorum 0,7 ; le calcul d'ecoPrimers est `floor(|E| x 0,7)` séquences (> longueur d'amorce). + +## 10. Mode A, étape 2b : mots stricts de bout en bout (2026-11-03) + +Code : `pkg/obiprimer/roles.go` (`Role`, `Selection`, `BelongingTo`) et `words.go` (`WordSearch`, +`StrictWords`, `FindStrictWords`, `MinSequences`). Tests : `words_test.go`. + +- `Selection` reproduit `-r` (restrict), `-i` (ignore) et `-E` (exception) d'ecoPrimers, vérifié sur les + quatre cas du tableau de la section 4 : exemple = dans Restrict et hors Except ; sinon contre-exemple si + dans Ignore ou Except ; sinon ignoré. Prédicats construits avec `obiseq.IsSubCladeOf` (existant). +- `FindStrictWords` lit les séquences une fois (consommateurs parallèles via `Split()`), indexe les + exemples dans un `SeqKmerCounter`, puis `CountAtLeast(floor(exemples x quorum))`. Quorum calculé avec + une petite tolérance (10 x 0,7 = 7) au lieu de la multiplication en float32 du C. + +**Contrôle contre ecoPrimers** (test opt-in : `OBIPRIMERS_TUTORIAL= go test ./pkg/obiprimer/`) : + +| Jeu | Exemples | Quorum | Mots obiprimers | Mots ecoPrimers | +|---|---|---|---|---| +| Teleostei | 3 909 | 2 736 | **255** | 255 (identique) | +| Cyprinidae | 325 | 227 | **1 207** | 1 208 | + +Écart de Cyprinidae : un mot en plus chez ecoPrimers, `TTTATGAAAACAAATAAA`. Comparaison des deux listes +complètes (le journal d'ecoPrimers imprime ses 1 208 mots) : tous les autres mots sont identiques. Un +comptage par force brute (deux brins, jonction circulaire) trouve ce mot dans **225** séquences sur 325, +sous le quorum de 227 : ecoPrimers le garde à tort. Le mécanisme exact n'est pas établi (`merge.c` élague +avec `remainingSeq + count >= quorum` pendant une fusion en place avec file d'attente ; non investigué +plus loin). Impact : ce mot n'apparaît dans **aucune** des 3 448 paires finales de +`Cyprinidae.ecoprimers`. Le test fige 1 207 avec cette explication. + +Mesure : étape « mots stricts » seule, 0,5 s (Cyprinidae) et 1,1 s (Teleostei) sur cette machine, lecture +de la taxonomie et du FASTA non comprise ; ecoPrimers complet : 134 s et 350 s (non décomposé par étape). + +**Bug existant repéré hors périmètre** : `go test -race` signale une course dans +`obiutils/pipe.go` (`globalLockerCounter` est un `int` incrémenté et décrémenté depuis plusieurs +goroutines, utilisé seulement pour des `log.Debugln`). Sans rapport avec obiprimers ; un `atomic.Int64` +suffirait. Non modifié (code partagé) : à décider. + +Suite : étape 3 (recherche approchée des mots stricts, <= e substitutions, deux brins, dans exemples +et contre-exemples), puis paires et indices. diff --git a/pkg/obiprimer/roles.go b/pkg/obiprimer/roles.go new file mode 100644 index 0000000..607399b --- /dev/null +++ b/pkg/obiprimer/roles.go @@ -0,0 +1,76 @@ +// Package obiprimer implements the core of obiprimers, the OBITools4 +// re-implementation of ecoPrimers: finding conserved words usable as primers +// in a set of example sequences, and the primer pairs they form. +package obiprimer + +import ( + "fmt" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitax" +) + +// Role tells how a sequence takes part in the primer search. +type Role uint8 + +const ( + // Ignored sequences are not used at all. + Ignored Role = iota + // Example sequences are those the primers must amplify. + Example + // CounterExample sequences are those the primers must not amplify. + CounterExample +) + +// Selection defines the roles from the taxonomy, with the ecoPrimers +// semantics of -r (restrict), -i (ignore) and -E (exception): +// +// - a sequence is an Example when it belongs to Restrict (all sequences when +// Restrict is nil) and not to Except; +// - otherwise it is a CounterExample when it belongs to Ignore or to Except; +// - otherwise it is Ignored. +// +// Nil predicates mean "no such constraint". +type Selection struct { + Restrict obiseq.SequencePredicate + Ignore obiseq.SequencePredicate + Except obiseq.SequencePredicate +} + +// Role returns the role of seq. +func (s Selection) Role(seq *obiseq.BioSequence) Role { + inRestrict := s.Restrict == nil || s.Restrict(seq) + inExcept := s.Except != nil && s.Except(seq) + + switch { + case inRestrict && !inExcept: + return Example + case inExcept, s.Ignore != nil && s.Ignore(seq): + return CounterExample + default: + return Ignored + } +} + +// BelongingTo returns a predicate true for the sequences whose taxon is a +// sub-clade of at least one of the given taxids. It returns nil (no +// constraint) when taxids is empty. +func BelongingTo(taxonomy *obitax.Taxonomy, taxids []string) (obiseq.SequencePredicate, error) { + var pred obiseq.SequencePredicate + for _, id := range taxids { + taxon, _, err := taxonomy.Taxon(id) + if err != nil { + return nil, fmt.Errorf("obiprimer: taxid %q: %w", id, err) + } + if taxon == nil { + return nil, fmt.Errorf("obiprimer: taxid %q is unknown from the taxonomy", id) + } + p := obiseq.IsSubCladeOf(taxonomy, taxon) + if pred == nil { + pred = p + } else { + pred = pred.Or(p) + } + } + return pred, nil +} diff --git a/pkg/obiprimer/words.go b/pkg/obiprimer/words.go new file mode 100644 index 0000000..c154e4d --- /dev/null +++ b/pkg/obiprimer/words.go @@ -0,0 +1,118 @@ +package obiprimer + +import ( + "fmt" + "math" + "sync" + "sync/atomic" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiiter" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer" +) + +// WordSearch holds the parameters of the strict word search. +type WordSearch struct { + PrimerLength int // word size (ecoPrimers -O, 18 by default; 2..31) + Quorum float64 // fraction of example sequences that must contain the word (ecoPrimers -q, 0.7) + Circular bool // sequences are circular (ecoPrimers -c) + Partitions int // partition files of the counter (default 256) + TmpDir string // parent of the temporary directory ("" = system default) +} + +// Counts gives the number of sequences by role. +type Counts struct { + Examples int + CounterExamples int + Ignored int +} + +// StrictWords is the result of the strict word search: the canonical words of +// PrimerLength bases present, without any error, in at least MinSequences +// example sequences. +type StrictWords struct { + PrimerLength int + MinSequences int // quorum in number of example sequences + Words []obikmer.SeqKmerCount + Counts Counts +} + +// MinSequences returns the quorum in sequences: floor(examples * quorum). +// The tiny epsilon makes exact products (10 * 0.7 = 7) robust to the binary +// representation of the fraction. +func MinSequences(examples int, quorum float64) int { + return int(math.Floor(float64(examples)*quorum + 1e-9)) +} + +// FindStrictWords reads every sequence of seqs once, indexes the example +// sequences, and returns the words reaching the quorum. Sequences no longer +// than the word size are counted as examples but cannot contain any word, as +// in ecoPrimers. Words are canonical (the smaller of the word and its reverse +// complement), so both strands are searched. +func FindStrictWords(seqs obiiter.IBioSequence, sel Selection, cfg WordSearch) (*StrictWords, error) { + if cfg.PrimerLength < 2 || cfg.PrimerLength > 31 { + return nil, fmt.Errorf("obiprimer: primer length must be between 2 and 31, got %d", cfg.PrimerLength) + } + if cfg.Quorum <= 0 || cfg.Quorum > 1 { + return nil, fmt.Errorf("obiprimer: quorum must be in ]0,1], got %g", cfg.Quorum) + } + if cfg.Partitions <= 0 { + cfg.Partitions = 256 + } + + counter, err := obikmer.NewSeqKmerCounter(cfg.TmpDir, cfg.PrimerLength, 0, cfg.Partitions) + if err != nil { + return nil, err + } + defer counter.Close() + + var nextID atomic.Uint32 + var examples, counters, ignored atomic.Int64 + + consumer := func(it obiiter.IBioSequence, wg *sync.WaitGroup) { + defer wg.Done() + for it.Next() { + for _, seq := range it.Get().Slice() { + switch sel.Role(seq) { + case Example: + examples.Add(1) + if seq.Len() > cfg.PrimerLength { + counter.AddSequence(nextID.Add(1)-1, seq.Sequence(), cfg.Circular) + } + case CounterExample: + counters.Add(1) + default: + ignored.Add(1) + } + } + } + } + + var wg sync.WaitGroup + for i := 1; i < obidefault.ParallelWorkers(); i++ { + wg.Add(1) + go consumer(seqs.Split(), &wg) + } + wg.Add(1) + go consumer(seqs, &wg) + wg.Wait() + + res := &StrictWords{ + PrimerLength: cfg.PrimerLength, + Counts: Counts{ + Examples: int(examples.Load()), + CounterExamples: int(counters.Load()), + Ignored: int(ignored.Load()), + }, + } + res.MinSequences = MinSequences(res.Counts.Examples, cfg.Quorum) + if res.Counts.Examples == 0 { + return nil, fmt.Errorf("obiprimer: no example sequence selected") + } + + res.Words, err = counter.CountAtLeast(max(res.MinSequences, 1)) + if err != nil { + return nil, err + } + return res, nil +} diff --git a/pkg/obiprimer/words_test.go b/pkg/obiprimer/words_test.go new file mode 100644 index 0000000..1fbdb0d --- /dev/null +++ b/pkg/obiprimer/words_test.go @@ -0,0 +1,223 @@ +package obiprimer + +import ( + "math/rand" + "os" + "path/filepath" + "strings" + "testing" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiformats" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiiter" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" +) + +func idPrefix(prefix string) obiseq.SequencePredicate { + return func(s *obiseq.BioSequence) bool { return strings.HasPrefix(s.Id(), prefix) } +} + +func TestSelectionRoles(t *testing.T) { + seq := func(id string) *obiseq.BioSequence { return obiseq.NewBioSequence(id, []byte("acgt"), "") } + // ids: "r" restricted clade, "x" exception inside it, "i" ignored clade, "o" other + cases := []struct { + name string + sel Selection + want map[string]Role + }{ + {"nothing set: everything is an example", + Selection{}, map[string]Role{"r1": Example, "x1": Example, "i1": Example, "o1": Example}}, + {"restrict only: the rest is dropped", + Selection{Restrict: idPrefix("r")}, + map[string]Role{"r1": Example, "x1": Ignored, "i1": Ignored, "o1": Ignored}}, + {"restrict + ignore", + Selection{Restrict: idPrefix("r"), Ignore: idPrefix("i")}, + map[string]Role{"r1": Example, "i1": CounterExample, "o1": Ignored}}, + {"exception inside restrict becomes a counterexample", + Selection{Restrict: idPrefix("r").Or(idPrefix("x")), Except: idPrefix("x")}, + map[string]Role{"r1": Example, "x1": CounterExample, "o1": Ignored}}, + {"example wins over ignore when both match", + Selection{Restrict: idPrefix("r"), Ignore: idPrefix("r").Or(idPrefix("i"))}, + map[string]Role{"r1": Example, "i1": CounterExample, "o1": Ignored}}, + } + for _, tc := range cases { + for id, want := range tc.want { + if got := tc.sel.Role(seq(id)); got != want { + t.Errorf("%s: %s: role %d, want %d", tc.name, id, got, want) + } + } + } +} + +func TestMinSequences(t *testing.T) { + for _, tc := range []struct { + n int + q float64 + want int + }{{10, 0.7, 7}, {3909, 0.7, 2736}, {325, 0.7, 227}, {100, 0.9, 90}, {7, 0.9, 6}, {1, 1, 1}} { + if got := MinSequences(tc.n, tc.q); got != tc.want { + t.Errorf("MinSequences(%d,%g) = %d, want %d", tc.n, tc.q, got, tc.want) + } + } +} + +func revcomp(s string) string { + c := map[byte]byte{'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'} + out := make([]byte, len(s)) + for i := range s { + out[len(s)-1-i] = c[s[i]] + } + return string(out) +} + +func randDNA(rng *rand.Rand, n int) string { + b := make([]byte, n) + for i := range b { + b[i] = "ACGT"[rng.Intn(4)] + } + return string(b) +} + +func TestFindStrictWords(t *testing.T) { + const k = 18 + rng := rand.New(rand.NewSource(1)) + motif := randDNA(rng, 30) + + var data obiseq.BioSequenceSlice + add := func(id, s string) { data = append(data, obiseq.NewBioSequence(id, []byte(s), "")) } + + // 10 examples; the motif is in 8 of them, 3 on the reverse strand. One + // example carries it twice. Another has an ambiguity breaking it. + for i := 0; i < 8; i++ { + m := motif + if i%3 == 0 { + m = revcomp(motif) + } + s := randDNA(rng, 120) + m + randDNA(rng, 120) + if i == 5 { + s += randDNA(rng, 50) + motif + } + add("ex"+string(rune('a'+i)), s) + } + broken := motif[:12] + "N" + motif[13:] + add("exh", randDNA(rng, 100)+broken+randDNA(rng, 100)) + add("exi", randDNA(rng, 300)) + // counterexamples and ignored sequences must not count + add("cx1", randDNA(rng, 50)+motif) + add("zz1", randDNA(rng, 50)+motif) + + sel := Selection{Restrict: idPrefix("ex").Or(idPrefix("cx")), Ignore: idPrefix("cx"), Except: idPrefix("cx")} + res, err := FindStrictWords(obiiter.IBatchOver("t", data, 3), sel, + WordSearch{PrimerLength: k, Quorum: 0.7, Partitions: 8, TmpDir: t.TempDir()}) + if err != nil { + t.Fatal(err) + } + + if res.Counts != (Counts{Examples: 10, CounterExamples: 1, Ignored: 1}) { + t.Fatalf("counts = %+v", res.Counts) + } + if res.MinSequences != 7 { + t.Fatalf("quorum = %d sequences, want 7", res.MinSequences) + } + + want := map[uint64]bool{} + for i := 0; i+k <= len(motif); i++ { + want[obikmer.EncodeCanonicalKmer([]byte(motif[i:i+k]), k)] = true + } + // the N breaks the windows that cover it, but 8 sequences still carry all of them + if len(res.Words) != len(want) { + t.Fatalf("%d words, want %d (the motif windows)", len(res.Words), len(want)) + } + for _, w := range res.Words { + if !want[w.Kmer] { + t.Errorf("unexpected word %d", w.Kmer) + } + if w.Sequences != 8 { + t.Errorf("word %d in %d sequences, want 8", w.Kmer, w.Sequences) + } + } + repeated := 0 + for _, w := range res.Words { + if w.Repeated { + repeated++ + } + } + if repeated != len(want) { + t.Errorf("%d words flagged Repeated, want all %d (motif twice in one sequence)", repeated, len(want)) + } +} + +func TestFindStrictWordsParameters(t *testing.T) { + it := func() obiiter.IBioSequence { + return obiiter.IBatchOver("t", obiseq.BioSequenceSlice{obiseq.NewBioSequence("a", []byte("acgtacgtacgtacgtacgtacgt"), "")}, 1) + } + if _, err := FindStrictWords(it(), Selection{}, WordSearch{PrimerLength: 1, Quorum: 0.7}); err == nil { + t.Error("primer length 1 must fail") + } + if _, err := FindStrictWords(it(), Selection{}, WordSearch{PrimerLength: 18, Quorum: 0}); err == nil { + t.Error("quorum 0 must fail") + } + if _, err := FindStrictWords(it(), Selection{Restrict: idPrefix("zz")}, + WordSearch{PrimerLength: 18, Quorum: 0.7}); err == nil { + t.Error("no example must fail") + } +} + +// Opt-in end-to-end check against ecoPrimers on the "Designing fish primers" +// tutorial data. Set OBIPRIMERS_TUTORIAL to the directory holding +// ncbitaxo_20250211.tgz and mito.vert.fasta.gz. +// +// ecoPrimers (-e 3 -3 2 -l 30 -L 150 -c, k=18, strict quorum 0.7) reports +// "Strict primer count": 1208 for Cyprinidae (-r 7953, 325 examples) and 255 +// for Teleostei (-r 32443, 3909 examples). +// +// Teleostei matches exactly. For Cyprinidae we find 1207: ecoPrimers keeps +// TTTATGAAAACAAATAAA, which a brute-force scan (both strands, circular +// junction) finds in only 225 of the 325 examples, below its own quorum of +// 227. That word is in none of the 3448 pairs ecoPrimers finally retains, so +// the divergence has no effect on the reference results. +func TestStrictWordsAgainstEcoPrimers(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 string + taxid string + examples int + words int + }{ + {"Cyprinidae", "7953", 325, 1207}, // ecoPrimers: 1208, see above + {"Teleostei", "32443", 3909, 255}, + } { + t.Run(tc.name, func(t *testing.T) { + restrict, err := BelongingTo(taxonomy, []string{tc.taxid}) + if err != nil { + t.Fatal(err) + } + // The predicates read the taxon through the taxid annotation. + seqs, err := obiformats.ReadSequencesFromFile(filepath.Join(dir, "mito.vert.fasta.gz")) + if err != nil { + t.Fatal(err) + } + res, err := FindStrictWords(seqs, Selection{Restrict: restrict}, + WordSearch{PrimerLength: 18, Quorum: 0.7, Circular: true, TmpDir: t.TempDir()}) + if err != nil { + t.Fatal(err) + } + t.Logf("%+v min=%d words=%d", res.Counts, res.MinSequences, len(res.Words)) + if res.Counts.Examples != tc.examples { + t.Errorf("%d examples, ecoPrimers reports %d", res.Counts.Examples, tc.examples) + } + if len(res.Words) != tc.words { + t.Errorf("%d strict words, ecoPrimers reports %d", len(res.Words), tc.words) + } + }) + } +}