From d6dd39f878b5b56936adf17395acf84db62f50b5 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Tue, 3 Nov 2026 21:12:04 +1100 Subject: [PATCH] Introduce sequence roles and filtering logic Implement new data structures for corpus management, multi-match filtering, and barcode resolution. --- OBIPRIMERS_MEMORY_PROMPT.md | 149 ++++++++++ OBIPRIMERS_NOTES.md | 137 ++++++++- pkg/obiprimer/amplicons.go | 232 +++++++++++++++ pkg/obiprimer/corpus.go | 41 ++- pkg/obiprimer/match.go | 8 +- pkg/obiprimer/multimatch.go | 11 +- pkg/obiprimer/pairs.go | 3 + pkg/obiprimer/pairs_eval.go | 323 ++++++--------------- pkg/obiprimer/pairs_test.go | 35 ++- pkg/obiprimer/resolution.go | 133 +++++++++ pkg/obiprimer/resolution_test.go | 153 ++++++++++ pkg/obiprimer/roles.go | 13 +- pkg/obiprimer/search.go | 2 +- pkg/obiprimer/search_real_test.go | 55 ++-- pkg/obiprimer/search_test.go | 81 +++++- pkg/obiprimer/taxa.go | 5 +- pkg/obiprimer/words_test.go | 12 +- pkg/obitools/obiprimers/embedded_test.go | 2 +- pkg/obitools/obiprimers/obiprimers.go | 8 +- pkg/obitools/obiprimers/obiprimers_test.go | 6 +- pkg/obitools/obiprimers/options.go | 8 +- pkg/obitools/obiprimers/output.go | 8 +- 22 files changed, 1083 insertions(+), 342 deletions(-) create mode 100644 OBIPRIMERS_MEMORY_PROMPT.md create mode 100644 pkg/obiprimer/amplicons.go create mode 100644 pkg/obiprimer/resolution.go create mode 100644 pkg/obiprimer/resolution_test.go diff --git a/OBIPRIMERS_MEMORY_PROMPT.md b/OBIPRIMERS_MEMORY_PROMPT.md new file mode 100644 index 0000000..a9666dd --- /dev/null +++ b/OBIPRIMERS_MEMORY_PROMPT.md @@ -0,0 +1,149 @@ +# Prompt : optimisation de l'usage mémoire de `obiprimers` + +Tu travailles dans le dépôt Go `/Users/coissac/Sync/travail/__MOI__/GO/obitools4` (module +`git.metabarcoding.org/obitools/obitools4/obitools4`, Go : `/usr/local/go/bin/go`, non présent dans le PATH). +Ta mission : **réduire l'usage mémoire de la commande `obiprimers`** sans changer ses résultats, en procédant +par la mesure. + +## 1. Contexte + +`obiprimers` est la ré-implémentation OBITools4 d'ecoPrimers (design de paires d'amorces à partir de génomes +entiers, avec évaluation de la qualité du marqueur délimité). Le mode « voie A » est terminé et validé. +Avant toute chose, **lis `OBIPRIMERS_NOTES.md`** (racine du dépôt) : la section « État courant », puis les +sections 11 à 13 (appariement approché), 16 et 16 bis (taxonomie embarquée ; **une erreur de mesure à ne pas +refaire**), 25 (définitions de Bc, Bs, Bd et élagage de la phase 2). + +Code concerné : +- `cmd/obitools/obiprimers/main.go`, `pkg/obitools/obiprimers/` (options, lecture, sortie CSV) ; +- `pkg/obiprimer/` : `corpus.go` (`LoadCorpus`), `taxa.go`, `words.go` (`FindStrictWords`), `match.go` et + `match_seed.go` (`MatchWords`), `pairs.go` (`selectPrimers`, `countPairs`), `pairs_eval.go` (`FindPairs`), + `amplicons.go` (`collectAmplicons`), `resolution.go`, `search.go` (`Search`) ; +- `pkg/obikmer/seqquorum*.go` (`SeqKmerCounter`, comptage par partitions sur disque) ; +- `pkg/obithermo/` (Tm), code partagé : `pkg/obiformats`, `pkg/obiiter`, `pkg/obiseq`, `pkg/obitax`. + +## 2. Décisions de l'utilisateur à respecter (non négociables) + +- **Tout est chargé en mémoire** : toutes les séquences de la base (cibles, contre-exemples, non-cibles) sont + gardées, parce que la lecture sur l'**entrée standard** doit continuer à marcher comme pour tout obitool. + Tu peux **compacter** ce qui est gardé ; tu ne peux pas passer à un schéma qui exige de relire le fichier. +- **Un seul calcul par indice** : ne touche pas aux définitions de Bc, Bs et Bd (`resolution.go`). Pas de + colonne ou d'option alternative. +- **Les résultats ne changent pas.** Critère : la sortie CSV doit rester **identique octet pour octet** sur les + cas de référence du §4. Si tu penses qu'un changement de résultat est justifié, **demande avant**. +- Aucune commande git qui modifie quoi que ce soit (pas de commit, add, mv, rm, reset…). Lecture seule + (`status`, `diff`, `log`) autorisée. Déplace les fichiers avec `mv`/`Write`/`Edit`. +- Fichiers source de moins d'environ 500 lignes. Pas de `.collect()` d'une grande structure quand un flux + suffit. Réutilise l'existant avant d'écrire du code ; enrichis une brique partagée plutôt que d'écrire du + code ad hoc, **mais** toute modification d'un paquet partagé (obiformats, obiiter, obiseq, obitax) doit + être signalée clairement et validée par les tests de ces paquets. +- **Pas de spéculation présentée comme un fait** : ce qui n'est pas mesuré est une hypothèse, dis-le. + Une régression de justesse est pire qu'un gain de mémoire : avant d'échanger l'une contre l'autre, demande. + +## 3. Ce qui est mesuré aujourd'hui (macOS, M3 Max, 16 cœurs) + +Jeu : 7 822 génomes mitochondriaux de vertébrés (131 Mb de bases), `-r 7953` (Cyprinidae : 325 cibles, 7 497 +non-cibles). Commande : +`obiprimers -r 7953 -e 3 --three-prime 2 -l 30 -L 150 -c -o out.csv` → 3 247 paires +(3 248 lignes avec l'en-tête). + +| Mesure | Valeur | +|---|---| +| pic RSS, taxonomie embarquée, fichier | **1 217 à 1 372 Mo** (bruit d'environ ±100 Mo d'une exécution à l'autre), 4,7 à 5,3 s | +| pic RSS, taxonomie embarquée, entrée standard | médiane **1 334 à 1 346 Mo** sur 8 exécutions entrelacées | +| tas des 7 822 séquences chargées | **193,7 Mo** (dont `taxonomic_path` : 19,2 Mo, soit 2,46 Ko par séquence) | +| même exécution **sans** l'élagage de phase 2 | 5 604 Mo (c'est ce qu'on a corrigé, voir §25 des notes) | +| avec `-t` (taxonomie NCBI complète) | +1,07 Go pour le seul chargeur (mesuré avec `obitaxonomy`), pic ~2,8 Go (mesuré avant les nouvelles définitions) | +| `GOGC=25` (ancienne mesure, `-t` NCBI) | pic 2,82 Go → 1,98 Go, mais temps CPU utilisateur ×2 | +| ecoPrimers (C) sur les mêmes jeux | 1,86 Go (Cyprinidae), 2,11 Go (Teleostei) | + +Tailles d'enregistrements : `ampRecord` 24 octets, `Match` 12, `seqMatch` 12, `hit` (intermédiaire de +`matchCorpus`) 16. Environ 6,85 millions d'occurrences approchées quand toutes les séquences sont chargées +(cas `-i 1`). + +**Le constat qui motive la mission** : les séquences gardées pèsent ~194 Mo, le pic du processus est ~1,3 Go. +Les 1,1 Go restants ne sont **pas expliqués** (hypothèses ci-dessous, aucune n'est vérifiée). + +## 4. Garde-fous à mettre en place AVANT de modifier quoi que ce soit + +1. **Binaire de référence** : `cd /Users/coissac/Sync/travail/__MOI__/GO/obitools4 && + GOEXPERIMENT=simd /usr/local/go/bin/go build -o /tmp/obiprimers.baseline ./cmd/obitools/obiprimers` +2. **Données** (les fichiers de `/tmp` peuvent avoir disparu, regénère-les) : dossier du tutoriel + `D="$HOME/Documents/Enseignement/DNA metabarcoding spring school/ecoPrimers-Fish"` : `mito.vert.fasta.gz`, + `ncbitaxo_20250211.tgz`, et pour le test d'échelle `mito.all.fasta.gz` / `mito.one.fasta.gz`. Fichier annoté + (taxonomie embarquée, évite le chargeur NCBI) : + `/Users/coissac/travail/__MOI__/GO/obitools4/build/obiannotate -t "$D/ncbitaxo_20250211.tgz" + --taxonomic-path -Z -o /tmp/mito.vert.path.fasta.gz "$D/mito.vert.fasta.gz"` + (les binaires déjà compilés sont dans `/Users/coissac/travail/__MOI__/GO/obitools4/build/`, sans `Sync`). +3. **CSV de référence**, à produire avec le binaire de référence et à conserver : Cyprinidae `-r 7953`, + Teleostei `-r 32443`, un cas avec contre-exemples `-r 7953 -i 1`, un cas `-U`, un cas `-R NC_083213`, et la + lecture sur entrée standard (`gzcat fichier | obiprimers ...`). Après chaque modification : `cmp` octet pour + octet. +4. **Tests** : `go test ./pkg/obiprimer/ ./pkg/obitools/obiprimers/ ./pkg/obithermo/ ./pkg/obikmer/ -race + -count=1`, et les tests opt-in sur les données réelles (durent ~75 s hors `-race`) : + `OBIPRIMERS_TUTORIAL="$D" go test ./pkg/obiprimer/ -count=1`. +5. **Méthode de mesure** (le piège principal) : le pic RSS bouge de ±100 Mo entre deux exécutions identiques + (rythme du ramasse-miettes). Mesure **au moins 5 exécutions entrelacées avant/après et compare les + médianes** (`/usr/bin/time -l` : octets sur macOS). Pour attribuer la mémoire, ne te fie pas au RSS : utilise + des profils de tas (`--pprof` de la CLI sert `/debug/pprof/heap` sur `localhost:6060`, `net/http/pprof` est + importé ; ou `runtime/pprof.WriteHeapProfile` aux frontières d'étapes ; `go tool pprof -sample_index=inuse_space` + et `alloc_space` ; `GODEBUG=gctrace=1`). **Erreur déjà commise (voir la section 16 bis des notes)** : en + mesurant le tas libéré par la suppression d'un attribut, la tranche de séquences n'était plus référencée + après la boucle et le GC a libéré *tout* ; utilise `runtime.KeepAlive` sur ce que tu veux garder vivant. + +## 5. Pistes à évaluer (toutes des HYPOTHÈSES à confirmer par un profil, classées sans ordre) + +- **Retard du ramasse-miettes** (`GOGC=100` laisse le tas monter jusqu'à ~2× le vivant). Pistes : limite + mémoire douce (`debug.SetMemoryLimit`) ou `debug.SetGCPercent` ajustés par la commande après la lecture ; + mesurer le coût CPU (l'ancienne mesure avec `GOGC=25` doublait le temps utilisateur). +- **Lecture** (code partagé `obiformats`/`obiiter`) : 131 Mb sont analysés, avec une table d'attributs et des + chaînes JSON par séquence. L'option `--batch-mem` (lots limités en mémoire) existe déjà : l'évaluer. +- **Séquences gardées** : `BioSequence` garde 1 octet par base et une table d'attributs. Ne garder que + l'identifiant, la séquence et le `taxid` (retirer les autres attributs, comme déjà fait pour + `taxonomic_path`, dans `dropTaxonomicPath` de `pkg/obitools/obiprimers/obiprimers.go`) ; ou stocker la + séquence à 2 bits par base (~33 Mo pour 131 Mb) avec une liste d'exceptions pour les bases ambiguës. Ce + second point touche `forEachOccurrence`, `ampHash`, `windowAsIs`, `windowRevComp` et `MatchWords`, qui lisent + les octets : à ne tenter que si le profil le justifie. +- **Occurrences stockées deux fois** : `[]WordMatches` (`matchCorpus`, intermédiaire `hit` de 16 octets par + occurrence puis recopie) et `primerSet.perSeq` (`selectPrimers`). Vérifier à quel moment chacune est + libérée, construire directement les structures finales, présizer les tranches (les `append` en doublement + créent un pic transitoire). +- **Phase 2b** (`collectAmplicons` sur les seules paires survivantes) : `ampRecord` de 24 octets (16 possibles + en empaquetant le rôle), enregistrements par worker ; **traitement des paires survivantes par lots bornés** + si le nombre d'enregistrements estimé dépasse un budget (idée notée en 25.2). Cas dur : `-T 0` (aucune paire + élaguée), où les tests atteignent ~5,5 Go de tas. +- **`SeqKmerCounter`** (`pkg/obikmer/seqquorum*.go`) : environ 16 octets par occurrence de k-mer pour les + partitions en cours de traitement ; le nombre de partitions est **fixé à 256**, alors que le total de bases + est connu une fois le corpus chargé : le choisir automatiquement. +- **Lecture en deux passes** pour la taxonomie embarquée (`readSequencesAndTaxonomy`) : la première passe ne + garde que les chemins, mais analyse tout. +- **Chargeur de taxonomie NCBI** : 1,07 Go à lui seul ; c'est du code partagé (`obitax`), hors périmètre sauf + si tu y trouves une économie sûre ; sinon, documente-le dans l'aide (la taxonomie embarquée évite ce coût). + +## 6. Test d'échelle (jamais fait) + +Tout est validé sur des génomes mitochondriaux (5 à 131 Mb). Mesure la croissance sur un fichier plus gros +(`mito.all.fasta.gz`, 172 Mo compressés, avec `--restrict-to-taxon` sur un grand clade, p. ex. les Metazoa) : +temps et pic RSS en fonction de la taille, pour savoir ce qui grossit vite (corpus, occurrences, amplicons, +compteur de k-mers). + +## 7. Livrables + +1. Un **rapport de profil** : où va la mémoire à chaque étape (lecture, sélection, mots stricts, appariement, + paires, phase 2a/2b), avec les profils de tas et les chiffres, avant toute modification. +2. Les optimisations **classées par gain mesuré et par risque**, puis implémentées une par une, chacune avec + mesure avant/après (médianes d'au moins 5 exécutions entrelacées), `cmp` des CSV de référence, tests verts + sous `-race`. +3. Une **nouvelle section** à la fin de `OBIPRIMERS_NOTES.md` (ne réécris pas l'historique) : mesures, + décisions, pistes écartées et pourquoi, ce qui reste. Pas de commit. +4. Si une optimisation exige de changer un résultat, une décision ou un paquet partagé : **arrête-toi et + demande**. + +## 8. Pièges connus de l'environnement + +- zsh ne découpe pas une variable en arguments : utilise des tableaux (`A=(-r 7953 ...)` puis `"${A[@]}"`) + ou des scripts `bash`. +- macOS : `sed -i ''` (sans suffixe) ; `/usr/bin/time -l` donne le pic en **octets** (`maximum resident set size`). +- `OBIPRIMERS_WORKERS` n'est lue que par les tests opt-in ; la CLI utilise `--max-cpu`. +- Les avertissements C de `obiformats/fastseq_read.h` sont connus et sans gravité. +- `go test -race` sur les tests opt-in est lent (~90 s) ; réserve-le aux tests unitaires. +- Le test `obiutils/TestSetString` était instable : il est corrigé (tri des membres dans `Set.String`). diff --git a/OBIPRIMERS_NOTES.md b/OBIPRIMERS_NOTES.md index 81f35ec..bec957a 100644 --- a/OBIPRIMERS_NOTES.md +++ b/OBIPRIMERS_NOTES.md @@ -3,6 +3,32 @@ Ré-implémentation d'ecoPrimers dans obitools4. Ce fichier garde l'état des recherches, des mesures et de l'avancement. Dernière mise à jour : 2026-11-03. +## État courant (à lire en premier) + +**Les trois indices de qualité d'une paire** (définitions en vigueur depuis la section 25 ; elles remplacent +celles d'ecoPrimers, qui figurent encore dans les sections 1 à 18 comme historique) : + +| Indice | Définition | +|---|---| +| **Bc** | taxons cibles amplifiés / taxons cibles | +| **Bs** | taxons cibles **identifiés** / taxons cibles amplifiés. Un taxon cible est identifié si **tous** ses barcodes ne sont portés que par lui, parmi **toutes** les séquences amplifiées par la paire (cibles, non-cibles, contre-exemples). Un barcode partagé rend mal identifiés **les deux** taxons. | +| **Bd** (nouveau) | part des barcodes portés par les cibles qui mènent sans ambiguïté à **un seul** taxon | + +Un seul calcul pour chacun, aucune variante ni colonne alternative. Les trois sont dans **toutes** les sorties +CSV (`bc`, `bs`, `bd`). Toute séquence sans taxon au rang choisi est un taxon anonyme distinct (un barcode +qu'elle partage est ambigu). Le seuil `-T` reste à 0,6. + +**Rôles des séquences** : cible (exemple), contre-exemple, **non-cible** (toutes les autres, qui ne sont plus +écartées : elles sont chargées en mémoire et leurs amplicons comptent dans Bs et Bd, pas dans les quorums). +Pour exclure des séquences de l'analyse, filtrer l'entrée en amont (`obigrep`). + +**Conséquence pour les validations contre ecoPrimers (sections 12, 15, 17, 18)** : elles restent exactes sur +ce qu'elles comparent (paires, drapeaux good, séquences et taxons amplifiés, Bc, Yule, longueurs, GC, Tm), mais +**Bs et `identifiés` ne sont plus comparables** à ceux d'ecoPrimers ; ils sont validés autrement (section 25). + +Sections en attente de décision de l'utilisateur : voie B (21.4), statistique de Tm (22.3-a), nombre de partitions +automatique, mémoire, tests shell, documentation, test `-e 0` (ecoPrimers n'avait pas terminé). + Sources de référence : - Article : `~/Sync/travail/__MOI__/ecoprimers/Documents/Riaz et al. 2011 ... .pdf` (NAR 39:e145) - Code C : `~/Sync/travail/__MOI__/ecoprimers/src/` (binaire : `/usr/local/bin/ecoprimers`) @@ -1027,7 +1053,7 @@ qui le portent. Un taxon cible t est **bien identifié** s'il existe un amplicon 3. **Séquence non cible sans taxon au rang choisi** : la traiter comme un taxon anonyme distinct (prudent : elle empêche l'identification d'une cible qui partage son amplicon) ou l'ignorer ? Je propose le premier. -## 24. Règle d'identification de Bs, et l'indice BD (2026-11-03) +## 24. Règle d'identification de Bs, et l'indice Bd (2026-11-03) **Statut : spécification, aucun code.** Complète et corrige la section 23. @@ -1038,10 +1064,10 @@ qui le portent. Un taxon cible t est **bien identifié** s'il existe un amplicon 2. **Règle d'identification : une espèce est bien identifiée si TOUS ses barcodes ne sont associés qu'à elle.** Dès qu'un barcode est partagé entre deux espèces, **les deux** sont mal identifiées. Cette règle remplace celle de 23.4-2 (que j'avais proposé de garder). -3. **Nouvel indice BD (barcode discrimination)** proposé par l'utilisateur : la part des **barcodes** qui +3. **Nouvel indice Bd (barcode discrimination)** proposé par l'utilisateur : la part des **barcodes** qui mènent sans ambiguïté à une espèce cible. Exemple de l'utilisateur : A porte les barcodes A1 et A2, B porte A2 et B1. Ni A ni B ne sont bien identifiées (elles partagent A2), mais **A1 et B1 sont discriminants**, A2 - ne l'est pas. Ici Bs = 0/2 et BD = 2/3. (Avec la règle d'ecoPrimers on aurait Bs = 2/2.) + ne l'est pas. Ici Bs = 0/2 et Bd = 2/3. (Avec la règle d'ecoPrimers on aurait Bs = 2/2.) ### 24.2 ecoPrimers applique bien l'autre règle (vérifié) @@ -1051,7 +1077,7 @@ taxon est donc compté **dès qu'un** de ses amplicons lui est propre. Confirmat cette règle reproduit toutes les colonnes `identifiés` et Bs de plus de vingt tables d'ecoPrimers (sections 12, 15, 17, 18) ; (ii) recalcul indépendant avec `obipcr` + Python sur la paire n° 1 de `Cyprinidae.ecoprimers` (320 amplicons chez `obipcr`, 318 chez ecoPrimers qui impose `-3 2`) : règle « au moins un barcode unique » -239/312 = 0,766 contre 237/310 = 0,765 pour ecoPrimers ; règle « tous uniques » 238/312 = 0,763 ; **BD = +239/312 = 0,766 contre 237/310 = 0,765 pour ecoPrimers ; règle « tous uniques » 238/312 = 0,763 ; **Bd = 240/265 = 0,906** (265 barcodes distincts). **Portée de (ii)** : sur cette paire les deux règles ne diffèrent que d'**une** espèce sur 312, car les données du tutoriel ont un génome par espèce (4 espèces seulement portent plusieurs barcodes) : ce recalcul ne tranche pas entre les règles à lui seul ; la preuve est le code @@ -1064,7 +1090,7 @@ et les tables reproduites. - **Non-régression** : la règle changeant, l'égalité des colonnes `identifiés` et Bs avec ecoPrimers n'est plus atteignable, même avec un corpus restreint aux cibles et contre-exemples (annule la proposition de 23.3 sur ce point). Toutes les autres colonnes (séquences, taxons, Bc, Yule, longueurs, GC, Tm, drapeaux good) restent - comparables. Pour Bs, BD et `identifiés` : tests sur cas synthétiques (dont l'exemple A1/A2/B1 ci-dessus) et + comparables. Pour Bs, Bd et `identifiés` : tests sur cas synthétiques (dont l'exemple A1/A2/B1 ci-dessus) et **vérification indépendante** par `obipcr` + script, comme en 24.2, sur des paires choisies. - **Seuil `-T` (spécificité minimale, 0,6 par défaut)** : calibré sur l'ancienne règle, plus indulgente. Le défaut sera peut-être à revoir ; à décider avec la spécification complète. @@ -1074,9 +1100,108 @@ et les tables reproduites. ### 24.4 Questions restantes -1. **BD** : l'inclure dans la sortie, et dans le classement ? Dénominateur proposé : les barcodes distincts +1. **Bd** : l'inclure dans la sortie, et dans le classement ? Dénominateur proposé : les barcodes distincts portés par **au moins une cible** (un barcode propre à une non-cible ne « mène » pas à une espèce cible) ; comptage par barcode distinct, sans pondération par le nombre de séquences. 2. **Non-cible sans taxon au rang choisi** (23.4-3, proposé, non confirmé) : taxon anonyme distinct, donc barcode partagé = ambigu (prudent). 3. **`-T` par défaut** (voir ci-dessus). + +## 25. Mise en œuvre des nouvelles définitions de Bc, Bs et Bd (2026-11-03) + +Demande de l'utilisateur : mettre à jour le document et `obiprimers` pour que les trois indices soient calculés +selon leur nouvelle définition. Décisions prises dans la discussion (sections 23 et 24) : un seul Bs ; tout en +mémoire ; dénominateur de Bs = cibles amplifiées ; identification par **tous** les barcodes ; Bd ajouté et +présent dans toutes les sorties ; non-cible sans taxon = taxon anonyme (ambigu) ; seuil `-T` inchangé (0,6). + +### 25.1 Changements de code + +- **Rôles** (`roles.go`, `corpus.go`) : `Ignored` devient `Other` (non-cible) ; `LoadCorpus` conserve toutes les + séquences ; `Counts.Others` ; `Corpus.Ignored` supprimé. Les quorums (sensibilité, faux positifs) ne + comptent que cibles et contre-exemples (`SequencesMatched`, phase 1 de `countPairs` qui saute les + non-cibles) ; `Taxa.CounterExamples` ne compte que les contre-exemples. +- **`-U`** (`DropMultiMatch(wm, corpus)`) : ne considère que les cibles et contre-exemples, comme ecoPrimers + (la validation de la section 15 reste exacte) ; une non-cible qui contient deux sites ne fait pas écarter le mot. +- **`resolution.go`** (nouveau) : `resolve` calcule, à partir des amplicons d'une paire, les cibles + amplifiées, les contre-exemples amplifiés, les taxons identifiés, les barcodes portés par les cibles et ceux + qui sont discriminants. États de barcode : ambigu / anonyme / porté par une cible. `lenientIdentified` est + la borne supérieure décrite en 25.2. +- **`amplicons.go`** (nouveau) : enregistrement d'un amplicon par séquence amplifiée, **tout rôle confondu** + (hash 64 bits, paire, taxon, rôle ; 24 octets). Chaque worker trie ses enregistrements par paire ; ils ne + sont ni copiés ni regroupés (l'ancienne version les dupliquait) ; les Tm minimales sont des minima par + worker, fusionnés à la fin ; l'évaluation des paires est parallèle. +- **`pairs_eval.go`** : `PairStats.Bd` ; `evaluate` utilise `resolve`. +- **Sortie CSV** : colonne `bd` ajoutée juste après `bs`. L'aide de `-r` décrit les non-cibles ; le journal + indique « other sequences » et le nombre de paires qui passent la borne de Bs. + +### 25.2 Un problème de mémoire, et son remède exact + +Première version : le pic de mémoire passait de **1,36 Go à 5,6 Go** (Cyprinidae, `-r 7953`), car la phase 2 +gardait l'amplicon de toutes les séquences amplifiées pour les 25 699 paires qui passent les quorums, alors que +seules ~3 400 seront retenues. + +Remède, **sans aucun changement de résultat** : un taxon identifié par la nouvelle règle (tous ses barcodes +propres parmi tout ce qui est amplifié) l'est aussi par la règle d'ecoPrimers (au moins un barcode propre parmi les +seules cibles), avec le même nombre de taxons cibles amplifiés. Donc `Bs_nouveau <= Bs_borne`. La phase 2 est +scindée : **2a**, amplicons des seules cibles (comme avant : Tm minimales, sites de la référence, borne de +Bs) ; **2b**, amplicons de toutes les séquences **pour les seules paires dont la borne atteint le seuil**. +Propriété vérifiée par `TestLenientIdentifiedIsAnUpperBound` (3 000 tirages aléatoires avec taxons anonymes, +plusieurs rôles, enregistrements en plusieurs morceaux). + +Résultat (Cyprinidae `-r 7953 -e 3 --three-prime 2 -l 30 -L 150 -c`, M3 Max) : + +| | pic mémoire | temps | paires | +|---|---|---|---| +| avant (ancienne définition) | 1 363 Mo | 6,5 s | 3 448 | +| nouvelle définition, sans élagage | 5 604 Mo | non mesuré | 3 247 | +| **nouvelle définition, avec élagage** | **1 217 à 1 372 Mo** | **4,7 à 5,3 s** | **3 247** | + +42 506 paires candidates, 25 699 passent les quorums, **3 448** passent la borne (c'est exactement le nombre de +l'ancienne table), **3 247** passent le seuil avec la nouvelle règle. Mise en garde : avec `-T 0` aucune paire +n'est élaguée et la mémoire remonte (les tests à `-T 0` atteignent ~5,5 Go de tas, taxonomie NCBI comprise). +Piste non engagée si cela devient gênant : traiter les paires survivantes par lots bornés. + +### 25.3 Valeurs : ancienne contre nouvelle règle (Cyprinidae, paire n° 1) + +`CAAGCCAGCCACATAACC` / `TTCTATAACGGGTGATGC` : 310 taxons cibles amplifiés ; ecoPrimers : 237 identifiés, +Bs = 0,7645 ; obiprimers : **233** identifiés, **Bs = 0,7516**, **Bd = 0,8935**. L'écart est faible ici car la +base a un génome par espèce (peu de barcodes multiples) et que les non-cibles amplifiées par cette paire sont +peu nombreuses à partager un barcode ; il sera plus marqué sur des bases à plusieurs séquences par espèce. + +### 25.4 Validation indépendante (obipcr + Python, base entière, `--three-prime 0`, `-T 0`) + +Six paires (les deux premières de la sortie et quatre tirées au hasard), amplicons obtenus par `obipcr -e 3 +-l 31 -L 149 -c` sur les **7 822 séquences**, Bs et Bd recalculés par un script qui n'utilise aucun code d'obiprimers. +Les paires amplifient de 3 202 à 7 796 séquences de la base entière (dont des milliers de non-cibles) : + +| Paire | taxons cibles amplifiés | identifiés | Bs | Bd | +|---|---|---|---|---| +| 1 | 314 = 314 | 237 = 237 | 0,755 = 0,755 | 0,892 = 0,892 | +| 2 | 314 = 314 | 237 = 237 | 0,755 = 0,755 | 0,892 = 0,892 | +| 3 | 314 = 314 | 5 = 5 | 0,016 = 0,016 | 0,357 = 0,357 | +| 4 | 315 = 315 | 14 = 14 | 0,044 = 0,044 | 14/27 = 0,5185 = 0,5185 | +| 5 | 315 = 315 | 7 = 7 | 0,022 = 0,022 | 0,233 = 0,233 | +| 6 | 314 = 314 | 15 = 15 | 0,048 = 0,048 | 0,469 = 0,469 | + +(format « obipcr+Python = obiprimers ».) Accord exact sur les 6 paires. Ces paires ont des Bs très bas : avec +tout ce qui est amplifié, un amplicon est vite partagé avec une non-cible. C'est l'effet recherché de la nouvelle +définition. Limite : la paire de référence 1 reste un cas où les deux règles diffèrent peu ; le cas qui les +distingue nettement (plusieurs barcodes par espèce) est couvert par les tests synthétiques (exemple A1/A2/B1 de +l'utilisateur, non-cible et contre-exemple partageant un barcode, taxon anonyme). + +### 25.5 Tests + +Nouveaux : `resolution_test.go` (exemple de l'utilisateur : Bs = 0/2, Bd = 2/3 ; tous uniques ; même espèce +plusieurs séquences ; non-cible ou contre-exemple partageant un barcode ; comptes de taxons ; taxon anonyme ; +ordre et découpage en morceaux ; borne supérieure) ; `TestNonTargetSharingABarcodeLowersBsAndBd` (de bout en +bout : une non-cible copie l'amplicon d'une cible, Bs = 7/8 et Bd = 7/8, quorums inchangés ; puis deux +contre-exemples amplifiés : 6 identifiés, Bd = 0,75, taxons contre-exemples = 2). + +Adaptés : les tests opt-in contre ecoPrimers (`TestPairsAgainstEcoPrimers`, `TestNoMultiMatchAgainstEcoPrimers`, +`TestCounterExamplesAgainstEcoPrimers`) désactivent le seuil de spécificité, vérifient que **chaque paire +d'ecoPrimers est trouvée** avec les mêmes comptes de séquences et de taxons, Yule, Bc, longueurs, GC, Tm et +drapeaux good, et que les taxons identifiés ne dépassent jamais ceux de la règle indulgente d'ecoPrimers. Résultats : +Cyprinidae 3 448 paires, 0 absente, 0 champ différent ; Teleostei 407, 0 absente, 6 à une séquence de moins +(l'artefact du `N`, section 12) ; `-U` 1 366, 0 absente, 0 différente ; contre-exemples 173, 0 absente, 0 +différente. obiprimers produit en plus 22 251 (Cyprinidae) et 2 389 (Teleostei) paires qu'ecoPrimers avait écartées avec son +propre Bs. diff --git a/pkg/obiprimer/amplicons.go b/pkg/obiprimer/amplicons.go new file mode 100644 index 0000000..828bd8a --- /dev/null +++ b/pkg/obiprimer/amplicons.go @@ -0,0 +1,232 @@ +package obiprimer + +import ( + "cmp" + "math" + "slices" + "sync" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obithermo" +) + +// ampRecord is an amplicon of a retained pair, in one sequence. +type ampRecord struct { + hash uint64 // the barcode: see ampHash + pair int32 + taxon int32 // taxon of the sequence at the chosen rank, -1 when it has none + role Role +} + +// refHit is an amplicon of a retained pair on the reference sequence. +type refHit struct { + up, down seqMatch + swapped bool +} + +// 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 +} + +// amplicons holds what the second phase gathers on the retained pairs. +type amplicons struct { + parts [][]ampRecord // one slice per worker, each sorted by pair + minTm1 []float32 // by pair: lowest Tm of the first primer over the example sites + minTm2 []float32 + refs map[int32]refHit +} + +// pair returns the amplicons of one retained pair, as one slice per worker. +func (a *amplicons) pair(i int32) [][]ampRecord { + out := make([][]ampRecord, 0, len(a.parts)) + byPair := func(r ampRecord, p int32) int { return cmp.Compare(r.pair, p) } + for _, part := range a.parts { + lo, _ := slices.BinarySearchFunc(part, i, byPair) + hi, _ := slices.BinarySearchFunc(part, i+1, byPair) + if hi > lo { + out = append(out, part[lo:hi]) + } + } + return out +} + +// collectOptions selects what collectAmplicons gathers. +type collectOptions struct { + examplesOnly bool // skip the counterexamples and the other sequences + model *obithermo.Model // melting temperatures of the example sites (nil: none) + reference bool // record the amplicons on the reference sequence +} + +// collectAmplicons gathers, in parallel, one record per amplicon of a retained +// pair, in the sequences of the corpus (all of them, or the examples only), and +// the lowest melting temperatures over the example sites. The records are not +// copied nor regrouped: each worker keeps its own, sorted by pair. +func (ps *primerSet) collectAmplicons(c *Corpus, taxa *Taxa, p PairParams, index map[pairKey]int32, opt collectOptions) *amplicons { + model := opt.model + lp := p.PrimerLength + workers := obidefault.ParallelWorkers() + npairs := len(index) + + res := &licons{ + parts: make([][]ampRecord, workers), + minTm1: make([]float32, npairs), + minTm2: make([]float32, npairs), + refs: map[int32]refHit{}, // written by the single worker that reads the reference + } + inf := float32(math.Inf(1)) + for i := range npairs { + res.minTm1[i], res.minTm2[i] = inf, inf + } + var tmMu sync.Mutex + + // 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)) + } + + jobs := make(chan int) + var wg sync.WaitGroup + for w := range workers { + wg.Add(1) + go func() { + defer wg.Done() + site1, site2 := make([]byte, lp), make([]byte, lp) + lowTm1, lowTm2 := map[int32]float32{}, map[int32]float32{} + for s := range jobs { + raw := c.Seqs[s].Sequence() + n := len(raw) + role := c.Roles[s] + if opt.examplesOnly && role != Example { + continue + } + ps.forEachOccurrence(c, s, p, func(o occurrence) { + pi, ok := index[o.key] + if !ok { + return + } + if opt.reference && p.HasReference && s == p.Reference { + if old, seen := res.refs[pi]; !seen || o.up.pos < old.up.pos { + res.refs[pi] = refHit{o.up, o.down, o.swapped} + } + } + res.parts[w] = append(res.parts[w], ampRecord{ + hash: ampHash(raw, int(o.up.pos)+lp, int(o.down.pos), o.swapped), + pair: pi, + taxon: taxa.Of[s], + role: role, + }) + + if role != Example || model == nil { + return + } + fp := oriented[o.up.primer][b2i(!o.up.reverse)] + rp := oriented[o.down.primer][b2i(o.down.reverse)] + tmF, tmR := inf, inf + if windowAsIs(site1, raw, int(o.up.pos), n, p.Circular) { + if t, err := model.Tm(fp, site1); err == nil && !math.IsNaN(t) { + tmF = float32(t) + } + } + if windowRevComp(site2, raw, int(o.down.pos), n, p.Circular) { + if t, err := model.Tm(rp, site2); err == nil && !math.IsNaN(t) { + tmR = float32(t) + } + } + if o.swapped { + tmF, tmR = tmR, tmF // slot 1 is the downstream primer + } + lowTm1[pi] = min(lowOr(lowTm1, pi, inf), tmF) + lowTm2[pi] = min(lowOr(lowTm2, pi, inf), tmR) + }) + } + + slices.SortFunc(res.parts[w], func(a, b ampRecord) int { return cmp.Compare(a.pair, b.pair) }) + tmMu.Lock() + for pi, t := range lowTm1 { + res.minTm1[pi] = min(res.minTm1[pi], t) + } + for pi, t := range lowTm2 { + res.minTm2[pi] = min(res.minTm2[pi], t) + } + tmMu.Unlock() + }() + } + for s := range c.Seqs { + jobs <- s + } + close(jobs) + wg.Wait() + return res +} + +func lowOr(m map[int32]float32, k int32, def float32) float32 { + if v, ok := m[k]; ok { + return v + } + return def +} + +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 +} diff --git a/pkg/obiprimer/corpus.go b/pkg/obiprimer/corpus.go index 489bf36..c64ea34 100644 --- a/pkg/obiprimer/corpus.go +++ b/pkg/obiprimer/corpus.go @@ -15,29 +15,32 @@ import ( type Counts struct { Examples int CounterExamples int - Ignored int + Others int } -// Corpus is the set of sequences used by the primer search: the examples and -// the counterexamples, in input order. Sequences with the Ignored role are not -// kept (as in ecoPrimers), only counted. +// Corpus is the set of sequences of the primer search, in input order, each +// with its role: the examples (the targets), the counterexamples, and all the +// other sequences, which are loaded as well since the amplicons they give +// decide whether a target is identified. // // The primer search needs several passes over the sequences (words are // searched approximately in every one of them), so the corpus lives in memory. type Corpus struct { - Seqs []*obiseq.BioSequence - Roles []Role // Roles[i] is Example or CounterExample - Ignored int // number of sequences dropped by the selection + Seqs []*obiseq.BioSequence + Roles []Role // Roles[i] is the role of Seqs[i] } // Counts returns the number of sequences by role. func (c *Corpus) Counts() Counts { - n := Counts{Ignored: c.Ignored} + var n Counts for _, r := range c.Roles { - if r == Example { + switch r { + case Example: n.Examples++ - } else { + case CounterExample: n.CounterExamples++ + default: + n.Others++ } } return n @@ -49,31 +52,24 @@ type corpusEntry struct { role Role } -// LoadCorpus reads every sequence of seqs, once, and keeps those selected as -// examples or counterexamples. The order of the input is preserved whatever +// LoadCorpus reads every sequence of seqs, once, and keeps all of them with the +// role given by the selection. The order of the input is preserved whatever // the number of reader goroutines. func LoadCorpus(seqs obiiter.IBioSequence, sel Selection) (*Corpus, error) { var mu sync.Mutex var entries []corpusEntry - ignored := 0 consumer := func(it obiiter.IBioSequence, wg *sync.WaitGroup) { defer wg.Done() var local []corpusEntry - localIgnored := 0 for it.Next() { batch := it.Get() for i, seq := range batch.Slice() { - if role := sel.Role(seq); role == Ignored { - localIgnored++ - } else { - local = append(local, corpusEntry{batch.Order(), i, seq, role}) - } + local = append(local, corpusEntry{batch.Order(), i, seq, sel.Role(seq)}) } } mu.Lock() entries = append(entries, local...) - ignored += localIgnored mu.Unlock() } @@ -94,9 +90,8 @@ func LoadCorpus(seqs obiiter.IBioSequence, sel Selection) (*Corpus, error) { }) c := &Corpus{ - Seqs: make([]*obiseq.BioSequence, len(entries)), - Roles: make([]Role, len(entries)), - Ignored: ignored, + Seqs: make([]*obiseq.BioSequence, len(entries)), + Roles: make([]Role, len(entries)), } for i, e := range entries { c.Seqs[i], c.Roles[i] = e.seq, e.role diff --git a/pkg/obiprimer/match.go b/pkg/obiprimer/match.go index 9a14f72..8f4216b 100644 --- a/pkg/obiprimer/match.go +++ b/pkg/obiprimer/match.go @@ -138,7 +138,8 @@ func matchCorpus(c *Corpus, words []uint64, mk scannerMaker) []WordMatches { } // SequencesMatched returns, for the matches of one word, the number of -// distinct example and counterexample sequences it occurs in. +// distinct example and counterexample sequences it occurs in. The other +// sequences are not counted. func (w WordMatches) SequencesMatched(c *Corpus) (examples, counterExamples int) { last := ^uint32(0) for _, m := range w.Matches { @@ -146,9 +147,10 @@ func (w WordMatches) SequencesMatched(c *Corpus) (examples, counterExamples int) continue } last = m.Seq - if c.Roles[m.Seq] == Example { + switch c.Roles[m.Seq] { + case Example: examples++ - } else { + case CounterExample: counterExamples++ } } diff --git a/pkg/obiprimer/multimatch.go b/pkg/obiprimer/multimatch.go index c06713f..fdfce0a 100644 --- a/pkg/obiprimer/multimatch.go +++ b/pkg/obiprimer/multimatch.go @@ -18,15 +18,16 @@ func DropRepeatedWords(words []obikmer.SeqKmerCount) (kept []obikmer.SeqKmerCoun } // DropMultiMatch removes the words having two or more approximate occurrences, -// both strands counted, in at least one sequence of the corpus, whatever its -// role. It is the second half of the ecoPrimers -U option: the primers that -// can anneal at several places of a sequence are discarded. -func DropMultiMatch(wm []WordMatches) (kept []WordMatches, dropped int) { +// both strands counted, in at least one example or counterexample sequence. It +// is the second half of the ecoPrimers -U option: the primers that can anneal +// at several places of a sequence they must (or must not) amplify are +// discarded. The other sequences are not considered. +func DropMultiMatch(wm []WordMatches, c *Corpus) (kept []WordMatches, dropped int) { kept = make([]WordMatches, 0, len(wm)) for _, w := range wm { multi := false for i := 1; i < len(w.Matches); i++ { - if w.Matches[i].Seq == w.Matches[i-1].Seq { + if w.Matches[i].Seq == w.Matches[i-1].Seq && c.Roles[w.Matches[i].Seq] != Other { multi = true break } diff --git a/pkg/obiprimer/pairs.go b/pkg/obiprimer/pairs.go index eaa2142..8457285 100644 --- a/pkg/obiprimer/pairs.go +++ b/pkg/obiprimer/pairs.go @@ -221,6 +221,9 @@ func (ps *primerSet) countPairs(c *Corpus, p PairParams) map[pairKey]*pairAcc { defer wg.Done() acc := partial[w] for s := range jobs { + if c.Roles[s] == Other { + continue // only the quorums are counted here + } isExample := c.Roles[s] == Example ps.forEachOccurrence(c, s, p, func(o occurrence) { a := acc[o.key] diff --git a/pkg/obiprimer/pairs_eval.go b/pkg/obiprimer/pairs_eval.go index 51c5189..aee5b5d 100644 --- a/pkg/obiprimer/pairs_eval.go +++ b/pkg/obiprimer/pairs_eval.go @@ -21,11 +21,12 @@ type PairStats struct { 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 + ExampleTaxa int // target taxa amplified + CounterTaxa int // counterexample taxa amplified + Bc float32 // coverage: ExampleTaxa over all the target taxa + Identified int // target taxa all of whose barcodes belong to them alone (see resolution.go) Bs float32 // specificity: Identified over ExampleTaxa + Bd float32 // barcode discrimination: share of the target barcodes leading to a single taxon MinLen, MaxLen int MeanLen float64 Tm1, Tm2 float64 // melting temperature of the primers against their perfect target @@ -50,44 +51,7 @@ type PairInfo struct { GoodPrimers int // of which good ones CandidatePairs int // pairs delimiting an amplicon in at least one sequence RetainedPairs int // pairs passing the sensitivity and false positive quorums, before the Bs filter -} - -// 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 -} - -// refHit is an amplicon of a retained pair on the reference sequence. -type refHit struct { - up, down seqMatch - swapped bool -} - -// 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 + Candidates int // of which pairs whose Bs upper bound reaches the threshold (their amplicons are gathered in full) } // FindPairs searches the primer pairs from the approximate occurrences of the @@ -145,24 +109,63 @@ func FindPairs(c *Corpus, taxa *Taxa, wm []WordMatches, p PairParams, model *obi index[k] = int32(i) } - // Phase 2: the amplicons of the retained pairs only. - records, refs := ps.collectAmplicons(c, p, index, model) - byPair := make([][]ampRecord, len(kept)) - for _, r := range records { - byPair[r.pair] = append(byPair[r.pair], r) + // Phase 2a: the amplicons of the target sequences, which give the melting + // temperatures, the reference sites, and an upper bound of Bs. It is cheap + // and rejects most pairs on the specificity threshold. + first := ps.collectAmplicons(c, taxa, p, index, collectOptions{examplesOnly: true, model: model, reference: true}) + var survivors []int // indices in kept + for i := range kept { + targets, identified := lenientIdentified(first.pair(int32(i))...) + if targets > 0 && float32(identified)/float32(targets) >= float32(p.MinSpecificity) { + survivors = append(survivors, i) + } } + first.parts = nil + info.Candidates = len(survivors) + + // Phase 2b: for the survivors only, the amplicons of every sequence, whatever + // its role: they decide which target taxa are identified. + index2 := make(map[pairKey]int32, len(survivors)) + for j, i := range survivors { + index2[kept[i]] = int32(j) + } + amps := ps.collectAmplicons(c, taxa, p, index2, collectOptions{}) + + results := make([]*PairStats, len(survivors)) + jobs := make(chan int) + var wg sync.WaitGroup + for range obidefault.ParallelWorkers() { + wg.Add(1) + go func() { + defer wg.Done() + for j := range jobs { + i := survivors[j] + k := kept[i] + st, ok := ps.evaluate(taxa, counts, p, k, acc[k], resolve(amps.pair(int32(j))...), + first.minTm1[i], first.minTm2[i], model) + if !ok { + continue + } + if h, found := first.refs[int32(i)]; found { + // evaluate orders the primers by string; the amplicon is read from the first one + s1 := orientedWord(ps.primers[k.p1].Word, p.PrimerLength, k.a1) + s2 := orientedWord(ps.primers[k.p2].Word, p.PrimerLength, k.a2) + st.RefLocation, st.RefAmplicon = referenceSites(c.Seqs[p.Reference].Sequence(), h, p, s2 < s1) + } + results[j] = &st + } + }() + } + for j := range survivors { + jobs <- j + } + close(jobs) + wg.Wait() var out []PairStats - for i, k := range kept { - st, ok := ps.evaluate(c, taxa, p, k, acc[k], byPair[i], model) - if ok { - if h, found := refs[int32(i)]; found { - // evaluate orders the primers by string; the amplicon is read from the first one - s1 := orientedWord(ps.primers[k.p1].Word, p.PrimerLength, k.a1) - s2 := orientedWord(ps.primers[k.p2].Word, p.PrimerLength, k.a2) - st.RefLocation, st.RefAmplicon = referenceSites(c.Seqs[p.Reference].Sequence(), h, p, s2 < s1) - } - out = append(out, st) + for _, st := range results { + if st != nil { + out = append(out, *st) } } slices.SortFunc(out, func(a, b PairStats) int { @@ -181,174 +184,21 @@ func boolCmp(a, b bool) int { 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, map[int32]refHit) { - lp := p.PrimerLength - workers := obidefault.ParallelWorkers() - parts := make([][]ampRecord, workers) +// evaluate puts the indices of one retained pair in a PairStats. It reports +// false when the pair must be dropped (no target taxon amplified, or +// specificity under the threshold). +func (ps *primerSet) evaluate(taxa *Taxa, counts Counts, p PairParams, k pairKey, a *pairAcc, res resolution, + minTm1, minTm2 float32, model *obithermo.Model) (PairStats, bool) { - // 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()) - - // The reference sequence is handled by a single worker, which is the only - // writer of refs; the first (leftmost) amplicon of each pair is kept. - refs := map[int32]refHit{} - - 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 - } - if p.HasReference && s == p.Reference { - if old, seen := refs[pi]; !seen || o.up.pos < old.up.pos { - refs[pi] = refHit{o.up, o.down, o.swapped} - } - } - 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, refs -} - -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 { + if res.targetTaxa == 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), + ExampleTaxa: res.targetTaxa, + CounterTaxa: res.counterTaxa, + Identified: res.identified, MinLen: int(a.minLen), MaxLen: int(a.maxLen), MeanLen: float64(a.sumLen) / float64(a.amplicons), @@ -359,12 +209,21 @@ func (ps *primerSet) evaluate(c *Corpus, taxa *Taxa, p PairParams, k pairKey, a 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) + st.Bc = float32(res.targetTaxa) / float32(taxa.ExampleTaxa) + st.Bs = float32(res.identified) / float32(res.targetTaxa) + st.Bd = float32(res.discriminant) / float32(res.barcodes) if st.Bs < float32(p.MinSpecificity) { return PairStats{}, false } + return ps.finish(st, p, k, a, minTm1, minTm2, model), true +} + +// finish adds the primers, their composition and melting temperatures, and puts +// the two primers in string order. +func (ps *primerSet) finish(st PairStats, p PairParams, k pairKey, a *pairAcc, + minTm1, minTm2 float32, model *obithermo.Model) PairStats { + 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) @@ -375,23 +234,11 @@ func (ps *primerSet) evaluate(c *Corpus, taxa *Taxa, p PairParams, k pairKey, a 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(minTm1), 1) { + st.MinTm1 = float64(minTm1) } - if !math.IsInf(float64(m1), 1) { - st.MinTm1 = float64(m1) - } - if !math.IsInf(float64(m2), 1) { - st.MinTm2 = float64(m2) + if !math.IsInf(float64(minTm2), 1) { + st.MinTm2 = float64(minTm2) } } @@ -402,5 +249,5 @@ func (ps *primerSet) evaluate(c *Corpus, taxa *Taxa, p PairParams, k pairKey, a st.Tm1, st.Tm2 = st.Tm2, st.Tm1 st.MinTm1, st.MinTm2 = st.MinTm2, st.MinTm1 } - return st, true + return st } diff --git a/pkg/obiprimer/pairs_test.go b/pkg/obiprimer/pairs_test.go index 9160cd7..563b5a9 100644 --- a/pkg/obiprimer/pairs_test.go +++ b/pkg/obiprimer/pairs_test.go @@ -74,9 +74,18 @@ func readEcoRows(t *testing.T, path string) []ecoRow { // 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 +// Bs and the number of identified taxa are NOT compared: obiprimers defines +// them differently (a taxon is identified only if all its barcodes are its own, +// and every amplified sequence counts, targets or not), so they are lower. The +// specificity threshold is switched off so that the pairs ecoPrimers dropped +// with its own Bs are still produced, and only the pairs of ecoPrimers are +// checked. What is checked: that every pair of ecoPrimers is found, with the +// same good flags, sequence and taxon counts, Yule, Bc, lengths, GC, Tm and +// lowest Tm; and that identified taxa never exceed those of ecoPrimers' more +// lenient rule on the same amplified sequences. +// +// Cyprinidae is reproduced exactly on these columns for its 3448 pairs. +// 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 @@ -143,7 +152,7 @@ func TestPairsAgainstEcoPrimers(t *testing.T) { 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, + Sensitivity: 0.9, FalsePositive: 0.1, MinSpecificity: 0, Circular: true, }, obithermo.DefaultModel()) if err != nil { t.Fatal(err) @@ -186,7 +195,7 @@ func TestPairsAgainstEcoPrimers(t *testing.T) { // 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 + diverges := g.Examples != w.in || g.ExampleTaxa != w.intaxa if diverges { divergent++ } @@ -196,15 +205,15 @@ func TestPairsAgainstEcoPrimers(t *testing.T) { 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) + if !diverges && g.Identified > w.well { + note("identified", w, "got %d, more than the %d of ecoPrimers' lenient rule", 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 !near(float64(g.Yule), w.yule, tol) || !near(float64(g.Bc), w.bc, tol) { + note("indices", w, "yule/bc got %.4f/%.4f want %.3f/%.3f", g.Yule, g.Bc, w.yule, w.bc) } 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) @@ -220,17 +229,11 @@ func TestPairsAgainstEcoPrimers(t *testing.T) { 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", + t.Logf("missing from ours: %d; pairs of ours that ecoPrimers dropped with its Bs: %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 diff --git a/pkg/obiprimer/resolution.go b/pkg/obiprimer/resolution.go new file mode 100644 index 0000000..86d0c89 --- /dev/null +++ b/pkg/obiprimer/resolution.go @@ -0,0 +1,133 @@ +package obiprimer + +// The three quality indices of a pair, over the sequences it amplifies. +// +// Every sequence the pair amplifies gives a barcode (the amplicon between the +// primers); all of them are taken into account, whatever their role: examples +// (the targets), counterexamples and non-targets. +// +// - Bc, coverage: target taxa amplified / target taxa. +// - Bs, specificity: target taxa identified / target taxa amplified. A target +// taxon is identified when ALL its barcodes belong to it alone: as soon as a +// barcode is shared with another taxon, target or not, counterexample or +// not, both taxa are badly identified. +// - Bd, barcode discrimination: the share of the barcodes carried by the +// targets that lead without ambiguity to one taxon. +// +// A sequence with no taxon at the chosen rank is an anonymous taxon, distinct +// from every other: a barcode it shares with a taxon is ambiguous. + +// barcodeState is what is known of one barcode while the amplicons are scanned. +type barcodeState struct { + taxon int32 // 1 + the first taxon met, 0 while none + flags uint8 +} + +const ( + barcodeAmbiguous = 1 << iota // two distinct taxa carry it + barcodeAnonymous // a sequence with no taxon carries it + barcodeTarget // a target sequence with a taxon carries it +) + +// ambiguous reports whether the barcode does not lead to a single taxon. +func (b barcodeState) ambiguous() bool { + return b.flags&barcodeAmbiguous != 0 || (b.flags&barcodeAnonymous != 0 && b.taxon != 0) +} + +// resolution is the outcome of the identification of the taxa of one pair. +type resolution struct { + targetTaxa int // target taxa amplified + counterTaxa int // counterexample taxa amplified + identified int // target taxa whose barcodes all belong to them alone + barcodes int // distinct barcodes carried by the targets + discriminant int // of which lead to a single taxon +} + +// resolve computes the resolution of one pair from its amplicons. The slices +// need not be sorted; they are only read. +func resolve(parts ...[]ampRecord) resolution { + var res resolution + + state := map[uint64]barcodeState{} + targets := map[int32]struct{}{} + counters := map[int32]struct{}{} + + for _, part := range parts { + for _, r := range part { + st := state[r.hash] + switch { + case r.taxon < 0: + st.flags |= barcodeAnonymous + case st.taxon == 0: + st.taxon = r.taxon + 1 + case st.taxon != r.taxon+1: + st.flags |= barcodeAmbiguous + } + if r.taxon >= 0 { + switch r.role { + case Example: + st.flags |= barcodeTarget + targets[r.taxon] = struct{}{} + case CounterExample: + counters[r.taxon] = struct{}{} + } + } + state[r.hash] = st + } + } + + bad := map[int32]struct{}{} + for _, part := range parts { + for _, r := range part { + if r.role == Example && r.taxon >= 0 && state[r.hash].ambiguous() { + bad[r.taxon] = struct{}{} + } + } + } + + for _, st := range state { + if st.flags&barcodeTarget != 0 { + res.barcodes++ + if !st.ambiguous() { + res.discriminant++ + } + } + } + res.targetTaxa, res.counterTaxa = len(targets), len(counters) + res.identified = len(targets) - len(bad) + return res +} + +// lenientIdentified is an upper bound of resolution.identified, computed on +// the target sequences alone and with the lenient rule of ecoPrimers (a taxon +// counts as soon as ONE of its barcodes belongs to it alone, among the targets). +// +// A taxon identified by resolve has all its barcodes alone among every amplified +// sequence, hence among the targets: it has at least one such barcode and is +// counted here too. So identified <= the second result, with the same number of +// target taxa, which lets a pair be rejected on the specificity threshold without +// gathering the amplicons of the non-target sequences. +func lenientIdentified(parts ...[]ampRecord) (targetTaxa, identified int) { + owner := map[uint64]int32{} // barcode -> its taxon, or -1 when shared + targets := map[int32]struct{}{} + for _, part := range parts { + for _, r := range part { + if r.role != Example || r.taxon < 0 { + continue + } + targets[r.taxon] = struct{}{} + if prev, seen := owner[r.hash]; !seen { + owner[r.hash] = r.taxon + } else if prev != r.taxon { + owner[r.hash] = -1 + } + } + } + alone := map[int32]struct{}{} + for _, t := range owner { + if t >= 0 { + alone[t] = struct{}{} + } + } + return len(targets), len(alone) +} diff --git a/pkg/obiprimer/resolution_test.go b/pkg/obiprimer/resolution_test.go new file mode 100644 index 0000000..4f625c9 --- /dev/null +++ b/pkg/obiprimer/resolution_test.go @@ -0,0 +1,153 @@ +package obiprimer + +import ( + "math/rand" + "testing" +) + +// rec builds an amplicon record: barcode, taxon (-1: none at the rank), role. +func rec(barcode uint64, taxon int32, role Role) ampRecord { + return ampRecord{hash: barcode, taxon: taxon, role: role} +} + +func TestResolutionUserExample(t *testing.T) { + // A has the barcodes A1 and A2, B has A2 and B1. A2 is shared: neither A + // nor B is identified, but A1 and B1 lead to one species each. + const A, B = 0, 1 + const A1, A2, B1 = 11, 12, 21 + res := resolve([]ampRecord{ + rec(A1, A, Example), rec(A2, A, Example), + rec(A2, B, Example), rec(B1, B, Example), + }) + want := resolution{targetTaxa: 2, identified: 0, barcodes: 3, discriminant: 2} + if res != want { + t.Fatalf("got %+v, want %+v", res, want) + } +} + +func TestResolutionAllUnique(t *testing.T) { + res := resolve([]ampRecord{rec(1, 0, Example), rec(2, 1, Example), rec(3, 2, Example)}) + if res.identified != 3 || res.discriminant != 3 || res.barcodes != 3 || res.targetTaxa != 3 { + t.Errorf("all barcodes distinct: %+v", res) + } +} + +func TestResolutionSameTaxonSharesItsBarcodeFreely(t *testing.T) { + // several sequences of one species with the same barcode are not ambiguous, + // and several barcodes for one species are fine when none is shared + res := resolve([]ampRecord{ + rec(1, 0, Example), rec(1, 0, Example), rec(1, 0, Example), + rec(2, 0, Example), + }) + if res.identified != 1 || res.barcodes != 2 || res.discriminant != 2 { + t.Errorf("%+v", res) + } +} + +// The barcode of a target that is unique among the targets but shared with a +// non-target, or with a counterexample, does not identify it. +func TestResolutionNonTargetsCount(t *testing.T) { + for _, role := range []Role{Other, CounterExample} { + res := resolve([]ampRecord{ + rec(1, 0, Example), // unique among the targets... + rec(1, 5, role), // ...but shared with a non-target taxon + rec(2, 1, Example), // identified + rec(3, 6, role), // a barcode of non-targets only + }) + if res.targetTaxa != 2 || res.identified != 1 { + t.Errorf("role %d: %+v, want 2 target taxa and 1 identified", role, res) + } + // barcodes carried by the targets: 1 (ambiguous) and 2; barcode 3 leads to a non-target + if res.barcodes != 2 || res.discriminant != 1 { + t.Errorf("role %d: Bd counts %d/%d, want 1/2", role, res.discriminant, res.barcodes) + } + } +} + +func TestResolutionAmplifiedTaxaCounts(t *testing.T) { + res := resolve([]ampRecord{ + rec(1, 0, Example), rec(2, 1, CounterExample), rec(3, 1, CounterExample), + rec(4, 2, CounterExample), rec(5, 3, Other), + }) + if res.targetTaxa != 1 || res.counterTaxa != 2 { + t.Errorf("target taxa %d (want 1), counterexample taxa %d (want 2: others are not counted)", + res.targetTaxa, res.counterTaxa) + } +} + +// A sequence with no taxon at the rank is an anonymous taxon, distinct from all +// the others: a barcode shared with it is ambiguous. +func TestResolutionAnonymousTaxon(t *testing.T) { + res := resolve([]ampRecord{ + rec(1, 0, Example), rec(1, -1, Other), // shared with an unknown: ambiguous + rec(2, 1, Example), + }) + if res.identified != 1 || res.barcodes != 2 || res.discriminant != 1 { + t.Errorf("target sharing with an anonymous sequence: %+v", res) + } + // two anonymous sequences with the same barcode never count: no target carries it + res = resolve([]ampRecord{rec(7, -1, Other), rec(7, -1, Other), rec(8, 0, Example)}) + if res.barcodes != 1 || res.discriminant != 1 || res.identified != 1 { + t.Errorf("anonymous barcodes alone: %+v", res) + } + // a target with no taxon is not a target taxon, and is anonymous too + res = resolve([]ampRecord{rec(1, -1, Example), rec(1, 0, Example), rec(2, 0, Example)}) + if res.targetTaxa != 1 || res.identified != 0 { + t.Errorf("an anonymous target sequence sharing a barcode must make the taxon ambiguous: %+v", res) + } +} + +func TestResolutionSeveralPartsAndOrder(t *testing.T) { + // the amplicons come as one slice per worker, in any order + a := []ampRecord{rec(1, 0, Example), rec(2, 1, Example)} + b := []ampRecord{rec(2, 2, Other), rec(3, 1, Example)} + if res := resolve(a, b); res.targetTaxa != 2 || res.identified != 1 { + t.Errorf("taxon 1 shares barcode 2 with taxon 2: %+v", res) + } + if res := resolve(b, a); res.targetTaxa != 2 || res.identified != 1 { + t.Errorf("order must not matter: %+v", res) + } +} + +// ecoPrimers' lenient rule on the user's example: A and B each own one barcode, +// so both count (the new rule counts none), and the bound is what lets a pair +// be rejected early. +func TestLenientIdentifiedUserExample(t *testing.T) { + recs := []ampRecord{ + rec(11, 0, Example), rec(12, 0, Example), + rec(12, 1, Example), rec(21, 1, Example), + } + if targets, identified := lenientIdentified(recs); targets != 2 || identified != 2 { + t.Errorf("lenient rule: %d of %d, want 2 of 2", identified, targets) + } + if res := resolve(recs); res.identified != 0 { + t.Errorf("strict rule: %d identified, want 0", res.identified) + } +} + +// The pruning of phase 2 is only exact if the bound never falls under the +// true count, with the same number of target taxa. +func TestLenientIdentifiedIsAnUpperBound(t *testing.T) { + rng := rand.New(rand.NewSource(3)) + for trial := 0; trial < 3000; trial++ { + var recs []ampRecord + n := 1 + rng.Intn(40) + for range n { + role := []Role{Example, Example, Example, CounterExample, Other}[rng.Intn(5)] + recs = append(recs, rec(uint64(rng.Intn(12)), int32(rng.Intn(8))-1, role)) // taxon -1: none + } + // split in a few parts, as the workers do + cut := rng.Intn(len(recs) + 1) + res := resolve(recs[:cut], recs[cut:]) + targets, lenient := lenientIdentified(recs[:cut], recs[cut:]) + if targets != res.targetTaxa { + t.Fatalf("trial %d: %d target taxa for the bound, %d for resolve", trial, targets, res.targetTaxa) + } + if res.identified > lenient { + t.Fatalf("trial %d: identified %d exceeds the bound %d for %v", trial, res.identified, lenient, recs) + } + if res.discriminant > res.barcodes || res.identified > res.targetTaxa { + t.Fatalf("trial %d: inconsistent %+v", trial, res) + } + } +} diff --git a/pkg/obiprimer/roles.go b/pkg/obiprimer/roles.go index 607399b..2ae584e 100644 --- a/pkg/obiprimer/roles.go +++ b/pkg/obiprimer/roles.go @@ -14,9 +14,12 @@ import ( type Role uint8 const ( - // Ignored sequences are not used at all. - Ignored Role = iota - // Example sequences are those the primers must amplify. + // Other sequences are non-targets: nothing is required of them, so they + // count neither in the sensitivity nor in the false positive quorums, but + // the amplicons they give are taken into account to decide whether a target + // is identified (see Bs and Bd). + Other Role = iota + // Example sequences are the targets: the primers must amplify them. Example // CounterExample sequences are those the primers must not amplify. CounterExample @@ -28,7 +31,7 @@ const ( // - 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. +// - otherwise it is Other. // // Nil predicates mean "no such constraint". type Selection struct { @@ -48,7 +51,7 @@ func (s Selection) Role(seq *obiseq.BioSequence) Role { case inExcept, s.Ignore != nil && s.Ignore(seq): return CounterExample default: - return Ignored + return Other } } diff --git a/pkg/obiprimer/search.go b/pkg/obiprimer/search.go index 61d14fc..f4a5a8f 100644 --- a/pkg/obiprimer/search.go +++ b/pkg/obiprimer/search.go @@ -87,7 +87,7 @@ func Search(c *Corpus, taxa *Taxa, p SearchParams) ([]PairStats, *SearchInfo, er return nil, info, err } if p.NoMultiMatch { - wm, info.MultiMatch = DropMultiMatch(wm) + wm, info.MultiMatch = DropMultiMatch(wm, c) } for _, w := range wm { info.Matches += len(w.Matches) diff --git a/pkg/obiprimer/search_real_test.go b/pkg/obiprimer/search_real_test.go index 53b9466..3e047f4 100644 --- a/pkg/obiprimer/search_real_test.go +++ b/pkg/obiprimer/search_real_test.go @@ -64,10 +64,14 @@ func tutorialParams(t *testing.T) SearchParams { // ecoPrimers -U on the Cyprinidae tutorial case (-e 3 -3 2 -l 30 -L 150 -c -r 7953 -U) // keeps 1366 pairs. testdata/ecoprimers_Cyprinidae_U.pairs.tsv holds, for each, // the primers, the amplified example and counterexample sequences, the amplified -// example and counterexample taxa, and the identified taxa. +// example and counterexample taxa, and the identified taxa. The identified taxa +// are not compared (obiprimers defines them differently, see +// TestPairsAgainstEcoPrimers) and the specificity threshold is off, so that +// the pairs ecoPrimers dropped with its own Bs do not hide any of its pairs. func TestNoMultiMatchAgainstEcoPrimers(t *testing.T) { corpus, taxa, _ := loadTutorial(t) p := tutorialParams(t) + p.MinSpecificity = 0 p.NoMultiMatch = true pairs, info, err := Search(corpus, taxa, p) if err != nil { @@ -91,23 +95,28 @@ func TestNoMultiMatchAgainstEcoPrimers(t *testing.T) { want[p1+"/"+p2] = v } - if len(pairs) != len(want) { - t.Errorf("%d pairs, ecoPrimers -U reports %d", len(pairs), len(want)) - } - bad := 0 + ours := map[string]PairStats{} for _, g := range pairs { - v, ok := want[g.Primer1+"/"+g.Primer2] + ours[g.Primer1+"/"+g.Primer2] = g + } + bad, missing := 0, 0 + for key, v := range want { + g, ok := ours[key] if !ok { - t.Errorf("pair %s/%s is not in the ecoPrimers -U result", g.Primer1, g.Primer2) + if missing++; missing <= 3 { + t.Errorf("ecoPrimers -U pair %s is not found", key) + } continue } - got := [5]int{g.Examples, g.CounterExamples, g.ExampleTaxa, g.CounterTaxa, g.Identified} - if got != v { + got := [4]int{g.Examples, g.CounterExamples, g.ExampleTaxa, g.CounterTaxa} + if got != [4]int{v[0], v[1], v[2], v[3]} || g.Identified > v[4] { if bad++; bad <= 3 { - t.Errorf("pair %s/%s: got %v, ecoPrimers %v", g.Primer1, g.Primer2, got, v) + t.Errorf("pair %s: got %v (identified %d), ecoPrimers %v", key, got, g.Identified, v) } } } + t.Logf("%d pairs of ecoPrimers -U checked, %d missing, %d different; obiprimers has %d pairs (%d more, dropped by ecoPrimers' own Bs)", + len(want), missing, bad, len(pairs), len(pairs)-len(want)+missing) } var locRE = regexp.MustCompile(`(\d+)\.\.(\d+)`) @@ -195,7 +204,9 @@ func TestCounterExamplesAgainstEcoPrimers(t *testing.T) { if c := corpus.Counts(); c.Examples != 325 || c.CounterExamples != 7497 || taxa.CounterExamples != 7387 { t.Fatalf("counts %+v, %d counterexample species (ecoPrimers: 325, 7497, 7387)", c, taxa.CounterExamples) } - pairs, info, err := Search(corpus, taxa, tutorialParams(t)) + pc := tutorialParams(t) + pc.MinSpecificity = 0 // ecoPrimers' own Bs filter is not reproduced, see TestPairsAgainstEcoPrimers + pairs, info, err := Search(corpus, taxa, pc) if err != nil { t.Fatal(err) } @@ -222,22 +233,26 @@ func TestCounterExamplesAgainstEcoPrimers(t *testing.T) { want[p1+"/"+p2] = r } - if len(pairs) != len(want) { - t.Errorf("%d pairs, ecoPrimers reports %d", len(pairs), len(want)) - } - bad := 0 + ours := map[string]PairStats{} for _, g := range pairs { - w, ok := want[g.Primer1+"/"+g.Primer2] + ours[g.Primer1+"/"+g.Primer2] = g + } + bad, missing := 0, 0 + for key, w := range want { + g, ok := ours[key] if !ok { - t.Errorf("pair %s/%s is not in the ecoPrimers result", g.Primer1, g.Primer2) + if missing++; missing <= 3 { + t.Errorf("ecoPrimers pair %s is not found", key) + } continue } good := string("bG"[b2i(g.Good1)]) + string("bG"[b2i(g.Good2)]) - got := [5]int{g.Examples, g.CounterExamples, g.ExampleTaxa, g.CounterTaxa, g.Identified} - if got != w.v || good != w.good { + got := [4]int{g.Examples, g.CounterExamples, g.ExampleTaxa, g.CounterTaxa} + if got != [4]int{w.v[0], w.v[1], w.v[2], w.v[3]} || good != w.good || g.Identified > w.v[4] { if bad++; bad <= 3 { - t.Errorf("pair %s/%s: got %s %v, ecoPrimers %s %v", g.Primer1, g.Primer2, good, got, w.good, w.v) + t.Errorf("pair %s: got %s %v (identified %d), ecoPrimers %s %v", key, good, got, g.Identified, w.good, w.v) } } } + t.Logf("%d pairs of ecoPrimers checked, %d missing, %d different; obiprimers has %d pairs", len(want), missing, bad, len(pairs)) } diff --git a/pkg/obiprimer/search_test.go b/pkg/obiprimer/search_test.go index 9bcc92e..2324ef2 100644 --- a/pkg/obiprimer/search_test.go +++ b/pkg/obiprimer/search_test.go @@ -74,7 +74,7 @@ func TestSearchFindsThePlantedPair(t *testing.T) { if p == nil { t.Fatalf("the planted pair was not found among %d pairs (%+v)", len(pairs), info) } - if p.Examples != 8 || p.MinLen != 40 || p.MaxLen != 40 || p.Bs != 1 || p.Bc != 1 { + if p.Examples != 8 || p.MinLen != 40 || p.MaxLen != 40 || p.Bs != 1 || p.Bc != 1 || p.Bd != 1 { t.Errorf("unexpected indices %+v", *p) } if p.RefLocation != "" { @@ -183,9 +183,11 @@ func TestDropMultiMatch(t *testing.T) { {Word: 2, Matches: []Match{{0, 5, 0, false}, {0, 50, 1, false}}}, // two in sequence 0 {Word: 3, Matches: []Match{{2, 5, 0, false}, {2, 5, 0, true}}}, // both strands of one sequence {Word: 4}, // no match at all: kept + {Word: 5, Matches: []Match{{3, 1, 0, false}, {3, 8, 0, false}}}, // twice, but in a non-target: kept } - kept, dropped := DropMultiMatch(wm) - if dropped != 2 || len(kept) != 2 || kept[0].Word != 1 || kept[1].Word != 4 { + c := corpusOf([]Role{Example, Example, CounterExample, Other}, "a", "a", "a", "a") + kept, dropped := DropMultiMatch(wm, c) + if dropped != 2 || len(kept) != 3 || kept[0].Word != 1 || kept[1].Word != 4 || kept[2].Word != 5 { t.Errorf("kept %v, dropped %d", kept, dropped) } } @@ -197,3 +199,76 @@ func TestDropRepeatedWords(t *testing.T) { t.Errorf("kept %v, dropped %d", kept, dropped) } } + +// A non-target amplified with the very barcode of a target taxon makes this +// target badly identified, although its barcode is unique among the targets; +// it counts in neither the sensitivity nor the false positive quorums. +func TestNonTargetSharingABarcodeLowersBsAndBd(t *testing.T) { + s := newSynthetic(8, -1, 11) + + alone, _, err := Search(s.corpus, s.taxa, s.params()) + if err != nil { + t.Fatal(err) + } + p0 := find(alone, s.f, s.r) + if p0 == nil || p0.Bs != 1 || p0.Bd != 1 { + t.Fatalf("without the non-target: %+v", p0) + } + + // a copy of sequence 3, of another taxon, that is not a target + copy3 := append([]byte(nil), s.corpus.Seqs[3].Sequence()...) + s.corpus.Seqs = append(s.corpus.Seqs, obiseq.NewBioSequence("other", copy3, "")) + s.corpus.Roles = append(s.corpus.Roles, Other) + s.taxa.Of = append(s.taxa.Of, 100) + + pairs, _, err := Search(s.corpus, s.taxa, s.params()) + if err != nil { + t.Fatal(err) + } + p := find(pairs, s.f, s.r) + if p == nil { + t.Fatal("the planted pair must still be found") + } + if p.Examples != 8 || p.ExampleTaxa != 8 || p.Bc != 1 { + t.Errorf("the non-target must not change the quorums nor Bc: %+v", *p) + } + if p.Identified != 7 || p.Bs != 0.875 { + t.Errorf("target 3 shares its barcode with the non-target: identified %d, Bs %v (want 7, 0.875)", p.Identified, p.Bs) + } + if p.Bd != 0.875 { + t.Errorf("Bd %v, want 7 discriminant barcodes of 8", p.Bd) + } + if p.CounterTaxa != 0 || p.CounterExamples != 0 { + t.Errorf("a non-target is not a counterexample: %+v", *p) + } + + // Now counterexamples: two copies (of sequences 3 and 5) that the pair amplifies and + // one unrelated sequence that it does not. 2 of 3 amplified is under the 0.9 quorum. + s.corpus.Roles[len(s.corpus.Roles)-1] = CounterExample + copy5 := append([]byte(nil), s.corpus.Seqs[5].Sequence()...) + s.corpus.Seqs = append(s.corpus.Seqs, obiseq.NewBioSequence("counter5", copy5, "")) + s.corpus.Roles = append(s.corpus.Roles, CounterExample) + s.taxa.Of = append(s.taxa.Of, 101) + unrelated := []byte(strings.ToLower(randDNA(rand.New(rand.NewSource(5)), 300))) + s.corpus.Seqs = append(s.corpus.Seqs, obiseq.NewBioSequence("unrelated", unrelated, "")) + s.corpus.Roles = append(s.corpus.Roles, CounterExample) + s.taxa.Of = append(s.taxa.Of, 102) + + pq := s.params() + pq.FalsePositive = 0.9 + pairs, _, err = Search(s.corpus, s.taxa, pq) + if err != nil { + t.Fatal(err) + } + pc := find(pairs, s.f, s.r) + if pc == nil { + t.Fatal("the pair amplifies 2 of 3 counterexamples (< 0.9): it must be retained") + } + if pc.CounterExamples != 2 || pc.CounterTaxa != 2 { + t.Errorf("counterexamples amplified %d (taxa %d), want 2 and 2", pc.CounterExamples, pc.CounterTaxa) + } + if pc.Identified != 6 || pc.Bd != 0.75 { + t.Errorf("targets 3 and 5 share their barcode with a counterexample: identified %d (want 6), Bd %v (want 0.75)", + pc.Identified, pc.Bd) + } +} diff --git a/pkg/obiprimer/taxa.go b/pkg/obiprimer/taxa.go index bd3417c..89698b0 100644 --- a/pkg/obiprimer/taxa.go +++ b/pkg/obiprimer/taxa.go @@ -58,9 +58,10 @@ func (c *Corpus) AssignTaxa(taxonomy *obitax.Taxonomy, rank string) (*Taxa, erro t.nodes = append(t.nodes, taxon) } t.Of[i] = id - if c.Roles[i] == Example { + switch c.Roles[i] { + case Example: exSeen[id] = struct{}{} - } else { + case CounterExample: coSeen[id] = struct{}{} } } diff --git a/pkg/obiprimer/words_test.go b/pkg/obiprimer/words_test.go index 8abaa18..1d22835 100644 --- a/pkg/obiprimer/words_test.go +++ b/pkg/obiprimer/words_test.go @@ -27,18 +27,18 @@ func TestSelectionRoles(t *testing.T) { }{ {"nothing set: everything is an example", Selection{}, map[string]Role{"r1": Example, "x1": Example, "i1": Example, "o1": Example}}, - {"restrict only: the rest is dropped", + {"restrict only: the rest are non-targets", Selection{Restrict: idPrefix("r")}, - map[string]Role{"r1": Example, "x1": Ignored, "i1": Ignored, "o1": Ignored}}, + map[string]Role{"r1": Example, "x1": Other, "i1": Other, "o1": Other}}, {"restrict + ignore", Selection{Restrict: idPrefix("r"), Ignore: idPrefix("i")}, - map[string]Role{"r1": Example, "i1": CounterExample, "o1": Ignored}}, + map[string]Role{"r1": Example, "i1": CounterExample, "o1": Other}}, {"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}}, + map[string]Role{"r1": Example, "x1": CounterExample, "o1": Other}}, {"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}}, + map[string]Role{"r1": Example, "i1": CounterExample, "o1": Other}}, } for _, tc := range cases { for id, want := range tc.want { @@ -116,7 +116,7 @@ func TestFindStrictWords(t *testing.T) { t.Fatal(err) } - if res.Counts != (Counts{Examples: 10, CounterExamples: 1, Ignored: 1}) { + if res.Counts != (Counts{Examples: 10, CounterExamples: 1, Others: 1}) { t.Fatalf("counts = %+v", res.Counts) } if res.MinSequences != 7 { diff --git a/pkg/obitools/obiprimers/embedded_test.go b/pkg/obitools/obiprimers/embedded_test.go index 8f6f7ed..cef413e 100644 --- a/pkg/obitools/obiprimers/embedded_test.go +++ b/pkg/obitools/obiprimers/embedded_test.go @@ -59,7 +59,7 @@ func TestEmbeddedTaxonomyFromFiles(t *testing.T) { if err != nil { t.Fatal(err) } - if c := corpus.Counts(); c.Examples != 2 || c.Ignored != 1 { + if c := corpus.Counts(); c.Examples != 2 || c.Others != 1 { t.Errorf("restricting to family 5: %+v", c) } diff --git a/pkg/obitools/obiprimers/obiprimers.go b/pkg/obitools/obiprimers/obiprimers.go index 21d036e..59e6dc1 100644 --- a/pkg/obitools/obiprimers/obiprimers.go +++ b/pkg/obitools/obiprimers/obiprimers.go @@ -140,9 +140,9 @@ func CLIPrimerPairs(sequences obiiter.IBioSequence, taxonomy *obitax.Taxonomy) ( } counts := corpus.Counts() - log.Infof("%d example sequences (%d %s), %d counterexamples (%d %s), %d ignored", + log.Infof("%d example sequences (%d %s), %d counterexamples (%d %s), %d other sequences", counts.Examples, taxa.ExampleTaxa, CLIRank(), - counts.CounterExamples, taxa.CounterExamples, CLIRank(), counts.Ignored) + counts.CounterExamples, taxa.CounterExamples, CLIRank(), counts.Others) pairs, info, err := obiprimer.Search(corpus, taxa, params) if info != nil { @@ -153,8 +153,8 @@ func CLIPrimerPairs(sequences obiiter.IBioSequence, taxonomy *obitax.Taxonomy) ( } log.Infof("%d approximate occurrences; %d conserved primers, %d good", info.Matches, info.Pair.Primers, info.Pair.GoodPrimers) - log.Infof("%d candidate pairs, %d amplifying enough examples, %d retained", - info.Pair.CandidatePairs, info.Pair.RetainedPairs, len(pairs)) + log.Infof("%d candidate pairs, %d amplifying enough examples, %d possibly specific enough, %d retained", + info.Pair.CandidatePairs, info.Pair.RetainedPairs, info.Pair.Candidates, len(pairs)) } return pairs, err } diff --git a/pkg/obitools/obiprimers/obiprimers_test.go b/pkg/obitools/obiprimers/obiprimers_test.go index b31474f..e13cce4 100644 --- a/pkg/obitools/obiprimers/obiprimers_test.go +++ b/pkg/obitools/obiprimers/obiprimers_test.go @@ -12,7 +12,7 @@ func samplePairs() []obiprimer.PairStats { { Primer1: "ACACCGCCCGTCACTCTC", Primer2: "CCAAGTGCACCTTCCGGT", Good1: true, Good2: false, Examples: 3837, CounterExamples: 2, Yule: 0.98249996, ExampleTaxa: 3805, CounterTaxa: 1, - Bc: 0.98169, Identified: 2650, Bs: 0.69647, MinLen: 54, MaxLen: 140, MeanLen: 74.75123, + Bc: 0.98169, Identified: 2650, Bs: 0.69647, Bd: 0.90566, MinLen: 54, MaxLen: 140, MeanLen: 74.75123, Tm1: 62.46, Tm2: 60.74, MinTm1: 36.81, MinTm2: math.NaN(), GC1: 12, GC2: 11, }, {Primer1: "AAAA", Primer2: "TTTT", Tm1: math.NaN(), Tm2: math.NaN(), MinTm1: math.NaN(), MinTm2: math.NaN()}, @@ -42,8 +42,8 @@ func TestPairsToCSVIterator(t *testing.T) { if r["index"] != 1 || r["primer1"] != "ACACCGCCCGTCACTCTC" || r["good1"] != true || r["good2"] != false { t.Errorf("unexpected record %v", r) } - if r["yule"] != 0.9825 || r["bc"] != 0.9817 || r["bs"] != 0.6965 || r["mean_length"] != 74.75 { - t.Errorf("indices are not rounded as expected: %v %v %v %v", r["yule"], r["bc"], r["bs"], r["mean_length"]) + if r["yule"] != 0.9825 || r["bc"] != 0.9817 || r["bs"] != 0.6965 || r["bd"] != 0.9057 || r["mean_length"] != 74.75 { + t.Errorf("indices are not rounded as expected: %v %v %v %v %v", r["yule"], r["bc"], r["bs"], r["bd"], r["mean_length"]) } if r["tm1"] != 62.5 || r["tm2"] != 60.7 || r["tm1_min"] != 36.8 { t.Errorf("temperatures: %v %v %v", r["tm1"], r["tm2"], r["tm1_min"]) diff --git a/pkg/obitools/obiprimers/options.go b/pkg/obitools/obiprimers/options.go index 968ed22..13d37f4 100644 --- a/pkg/obitools/obiprimers/options.go +++ b/pkg/obitools/obiprimers/options.go @@ -51,9 +51,11 @@ func OptionSet(options *getoptions.GetOpt) { options.StringSliceVar(&_restrict, "restrict-to-taxon", 1, 1, options.Alias("r"), options.ArgName("TAXID"), - options.Description("Restricts the examples (the sequences the primers must amplify) to the "+ - "sequences belonging to this taxon. Other sequences are ignored, unless they are "+ - "counterexamples. Can be repeated. Without this option every sequence is an example.")) + options.Description("Restricts the examples (the targets, the sequences the primers must amplify) to "+ + "the sequences belonging to this taxon. The other sequences are non-targets: they are not "+ + "required to be amplified, but the barcodes they give are taken into account to decide "+ + "whether a target is identified (Bs, Bd). Can be repeated. Without this option every "+ + "sequence is an example.")) options.StringSliceVar(&_ignore, "ignore-taxon", 1, 1, options.Alias("i"), diff --git a/pkg/obitools/obiprimers/output.go b/pkg/obitools/obiprimers/output.go index c2f987d..af37b7a 100644 --- a/pkg/obitools/obiprimers/output.go +++ b/pkg/obitools/obiprimers/output.go @@ -18,9 +18,10 @@ var csvColumns = []string{ "examples", "counterexamples", // 11-12 sequences amplified "yule", // 13 "example_taxa", "counterexample_taxa", // 14-15 taxa amplified - "bc", // 16 coverage - "identified_taxa", // 17 - "bs", // 18 specificity + "bc", // 16 coverage: target taxa amplified / target taxa + "identified_taxa", // 17 target taxa all of whose barcodes belong to them alone + "bs", // 18 specificity: identified / amplified target taxa + "bd", // barcode discrimination (new): share of the target barcodes leading to one taxon "min_length", "max_length", "mean_length", // 19-21 } @@ -44,6 +45,7 @@ func pairRecord(index int, p obiprimer.PairStats) obiitercsv.CSVRecord { "yule": round(float64(p.Yule), 4), "example_taxa": p.ExampleTaxa, "counterexample_taxa": p.CounterTaxa, "bc": round(float64(p.Bc), 4), "identified_taxa": p.Identified, "bs": round(float64(p.Bs), 4), + "bd": round(float64(p.Bd), 4), "min_length": p.MinLen, "max_length": p.MaxLen, "mean_length": round(p.MeanLen, 2), } if p.RefLocation != "" {