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) + } + }) + } +}