diff --git a/OBIPRIMERS_NOTES.md b/OBIPRIMERS_NOTES.md index 0a9717e..81f35ec 100644 --- a/OBIPRIMERS_NOTES.md +++ b/OBIPRIMERS_NOTES.md @@ -586,3 +586,497 @@ course dans `obiiter.IBioSequence.RebatchBySize` (appelée par `CLIReadBioSequen réaffectait la variable capturée `iterator` (`iterator = iterator.SortBatches()`) pendant que la fonction appelante lisait `iterator.IsPaired()`. Corrigé avec une variable locale (`sorted`), comportement inchangé ; les tests de `obiiter`, `obiformats` et `obitools/...` passent. + +### 16 bis. Suppression de `taxonomic_path` après extraction (2026-11-03, demande de l'utilisateur) + +Une fois la taxonomie construite, l'attribut `taxonomic_path` (liste d'une vingtaine de chaînes par +séquence) ne sert plus : la recherche lit les taxons par le `taxid`. Il est supprimé dans les trois +chemins de lecture (`obiprimers.go`) : fichiers (worker `obiseq.AnnotatorToSeqWorker` sur la seconde +passe), entrée standard (en place sur la tranche gardée en mémoire, avant de la réémettre) et `-t` donné +avec un fichier annoté (les chemins sont inutiles). Réutilisation : `AnnotatorToSeqWorker` et +`BioSequence.DeleteAttribute` existants ; seul ajout partagé : la constante +`obiseq.TaxonomicPathAttribute`, utilisée par `SetPath`/`Path` à la place du littéral. Test : +`TestTaxonomicPathIsDroppedAfterExtraction` (le chemin disparaît, le `taxid` reste). + +**Mesures, sans gain visible sur le pic de mémoire du processus** : +- tas réellement libéré (7 822 séquences, 131 Mb de bases, `GC` forcé, tranche gardée vivante) : + **19,2 Mo sur 193,7 Mo (9,9 %), 2,46 Ko par séquence**. Un premier essai donnait « 99,9 % » : il était faux + (la tranche n'était plus référencée, donc toutes les séquences avaient été libérées) ; +- pic de mémoire du processus (entrée standard, 8 mesures entrelacées chacune) : médiane **1 334 Mo avant, + 1 346 Mo après**, dispersion > 250 Mo : aucune différence détectable. + +Interprétation : le gain est réel mais de l'ordre de 10 % des séquences gardées, soit ~1 % du pic ici ; il +deviendrait important sur de très grandes bases (à 2,46 Ko par séquence, ~2,5 Go pour un million de +séquences). **Le pic de 1,3 Go est sans commune mesure avec les 194 Mo de séquences gardées** : l'essentiel +du pic vient d'ailleurs (lecture/analyse des 131 Mb du FASTA, retard du ramasse-miettes : `GOGC=25` avait +fait passer le pic de 2,82 à 1,98 Go). C'est la piste à analyser pour la mémoire. + +## 17. Validation avec contre-exemples : `-r 7953 -i 1` (2026-11-03) + +Premier cas **avec contre-exemples** comparé à ecoPrimers (la section 1 de l'inventaire : jusque-là tous les +jeux de référence avaient 0 contre-exemple). Les 325 Cyprinidae sont les exemples (`-r 7953`) et toutes les +autres séquences de la base (7 497, 7 387 espèces) sont des contre-exemples (`-i 1`, la racine). Paramètres : +`-e 3 -3 2 -l 30 -L 150 -c`. (« -I 1 » demandé : ecoPrimers n'a pas d'option `-I` ; sa boucle `getopt` +ignore silencieusement une option inconnue, donc c'était `-i`.) + +| | ecoPrimers | obiprimers | +|---|---|---| +| exemples / contre-exemples | 325 (316 esp.) / 7 497 (7 387 esp.) | identique | +| quorum de contre-exemples | 749 / 7 497 | identique | +| paires retenues | **173** | **173** | +| temps | **2 728 s (45,5 min)**, 1 cœur | **1,9 s**, 16 workers | +| pic mémoire | 756 Mo | 781 Mo (taxonomie embarquée) | + +**Les deux tables sont identiques sur tous les champs** (173/173 paires, 0 différence : amorces, drapeaux +good, séquences et taxons amplifiés côté exemples ET contre-exemples, Yule, Bc, Bs, longueurs, GC, Tm). +Répartition des drapeaux good : 97 `bG` et 76 `Gb` (une seule amorce good par paire, comme attendu avec +des contre-exemples aussi nombreux) ; contre-exemples amplifiés par paire : de 264 à 513 séquences. +Notre journal : 6 850 691 occurrences approchées (contre 227 549 sans contre-exemples), 1 200 amorces +conservées dont **14 good**, 2 037 paires candidates, 200 amplifiant assez d'exemples, 173 retenues +(ecoPrimers : « Total pair count 2050 », écart de 13 attendu : il a un mot strict de plus, section 10). + +Fixé dans `testdata/ecoprimers_Cyprinidae_I1.pairs.tsv` ; test opt-in `TestCounterExamplesAgainstEcoPrimers`. +Cela couvre désormais : quorum de faux positifs, drapeau good avec contre-exemples, comptes de taxons +contre-exemples. Reste à comparer (inventaire section 1) : séquences linéaires, `--rank` autre qu'espèce, +`-e 0`, autres `-e` / `-O` / quorums, et `-E` (exception). + +## 18. Validation par variantes d'options, Cyprinidae `-r 7953` (2026-11-03) + +16 exécutions d'ecoPrimers lancées en parallèle (1 cœur chacune, 2 à 6 minutes sans contre-exemples), +comparées ligne à ligne (multi-ensembles de lignes, car des paires d'amorces AT répétées existent des deux +côtés) à `obiprimers` avec les options équivalentes. Scripts jetables : `/tmp/run_val*.sh`, +`/tmp/run_mine.sh`, `/tmp/compare_val.py` (non versionnés). + +**Identiques sur toutes les lignes et tous les champs** : + +| Variante | Lignes | +|---|---| +| séquences linéaires (sans `-c`) | 3 448 | +| `-e 1` | 2 031 | +| `-e 5` | 4 080 | +| `-O 20` | 2 288 | +| `-O 25` | 867 | +| `-q 0.8 -s 0.95` | 1 110 | +| `-q 0.5 -s 0.8` | 13 941 | +| `-3 0` | 5 405 | +| `-3 4` | 2 262 | +| `-l 10 -L 80` | 223 | +| `-E 7954` (Danio, 5 contre-exemples hors des exemples) | 11 | +| `-E 7954 -x 0.05 -s 0.8` | 16 | +| `-E 84644 -x 0.95` (Labeo : 304 exemples, 21 contre-exemples dans les exemples) | 54 | + +`-E 84644` seul et `-E 84644 -x 0.05 -s 0.8` donnent 0 paire des deux côtés : identiques mais sans valeur, +les Labeo étant amplifiés par toutes les amorces conservées ; d'où la variante `-x 0.95`. + +**`--rank genus` et `-T 0.3` : différents, mais entièrement expliqués** (par mesure) : +1. **118 et 119 paires seulement chez ecoPrimers** : toutes contiennent le mot strict en trop d'ecoPrimers + `TTTATGAAAACAAATAAA` (section 10 : présent dans 225 séquences, quorum 227). Aucune des paires absentes + ne l'ignore ; aucune ligne des variantes identiques ne le contient. Leurs comptes d'exemples vont de 294 + à 325 : l'hypothèse « effet de seuil à 293 » avait été testée et réfutée. +2. **52 et 65 lignes avec une séquence d'exemple de plus chez ecoPrimers** (exactement +1 exemple, +0 + contre-exemple à chaque fois) : l'artefact du `N` dans la vérification du 3' (section 12). +3. **9 écarts de Tm minimale** (mêmes comptes de séquences) : ecoPrimers décode les sites contenant une + base ambiguë avec le même registre à décalage ; nous ignorons ces sites (section 12). +Preuve : un binaire temporaire émulant exactement `extractSite` (3' et sites de Tm), sources restaurées +ensuite (diff identique aux sauvegardes), supprime **toutes** les différences de champs ; ne restent que +les 118/119 paires du mot en trop. + +## 19. Décisions de l'utilisateur et options `-b`, `-M`, `-k` (2026-11-03) + +- **Taxonomie embarquée : extraction stricte** (comme `obitag`) : comportement actuel conservé, rien à changer. +- **Test instable corrigé** : `obiutils.Set.String()` parcourait une map Go (ordre aléatoire) ; elle trie + désormais ses membres (entiers, flottants et chaînes par valeur, y compris types nommés ; repli sur la forme + texte pour les autres types). `Members()` est inchangé. Test `TestSetStringIsStable` ; `TestSetString` passe + 200 fois de suite (il échouait ~4 fois sur 6) et tout le paquet passe. +- **`-b`, `-M`, `-k`, `-p`** : lus dans `PrimerSets.c` (voir ci-dessous) ; porter ou abandonner : à trancher. + +### Ce que sont ces options dans ecoPrimers (lu dans le code) + +Un « lien » entre deux paires = le nombre de séquences que **les deux** paires identifient sans ambiguïté +(`wellIdentifiedSeqs[k] == 1` pour les deux). Le but est de choisir des paires complémentaires (plusieurs +marqueurs) : peu de liens = les paires identifient des taxons différents. +- `-b` : filtre la liste finale des paires avec `primers_changeSortedArray` (puis réimprime). Une paire est + gardée s'il existe une autre paire avec laquelle elle a au plus 5 liens ou plus de la moitié du plus petit + nombre de taxons identifiés ; elle n'est écartée que si, avec **toutes** les autres, ses liens sont entre 6 + et cette moitié. Règle heuristique sans justification dans le code. +- `-M` : seuil absolu de liens (nom trompeur : « percent » mais utilisé comme un nombre). Quand il est donné, + la boucle contient un `break` inconditionnel : chaque paire n'est comparée qu'à la **première** autre paire + de la liste, donc le filtre ne fait pas ce qu'il annonce. +- `-k` : n'est lu que par `ok_to_add` (construction d'ensembles d'amorces), inatteignable puisque `-p` est + désactivé : option morte. +- `-p` : sa ligne d'affectation est commentée dans `ecoprimer.c` ; elle ne fait rien. Le code d'ensembles + associé (glouton, recuit simulé, tabou, force brute) n'est jamais exécuté. +- `primers_filterWithGivenLinks` et `primers_graph_graphviz` ne sont jamais appelées. +Dépendance technique d'un éventuel portage : il faudrait garder, pour chaque paire retenue, l'ensemble des +séquences qu'elle identifie (aujourd'hui seule leur nombre de taxons est conservé). + +## 20. Le problème multi-marqueurs : spécification (discussion du 2026-11-03) + +**Statut : spécification en cours, aucun code écrit.** Décisions de l'utilisateur ci-dessous ; les points +marqués « proposition » ou « question ouverte » sont de moi et n'ont pas été validés. + +### 20.1 Origine de `-b`, `-M`, `-k`, `-p` : abandonnées + +Ces options sont des pistes explorées pendant la thèse de Tiayyba Riaz (auteur de l'article), il y a une +quinzaine d'années, jamais menées à terme : le problème était mal posé. L'utilisateur avait oublié qu'elles +traînaient dans le code. **Décision : elles ne sont pas portées** (confirmant l'analyse de la section 19 : deux +options mortes, une buguée, un filtre heuristique non justifié). + +Le modèle de la thèse : chaque paire d'amorces couvre et reconnaît une partie des espèces, avec un +recouvrement partiel entre paires ; on cherche la meilleure combinaison (graphe de « liens », glouton, recuit +simulé, tabou, force brute). Ce modèle est quasi inapplicable : un recouvrement d'espèces entre deux paires +visant des clades différents signifie qu'une amorce est partagée entre ces clades, donc que le recouvrement +est presque toujours **complet**. Les ensembles de couverture sont emboîtés ou disjoints par clade. + +### 20.2 Formulation correcte + +But : un ensemble de marqueurs (un multiplex) qui, ensemble, discriminent un maximum d'espèces. +Structure : +1. **une paire de base** qui couvre la totalité du groupe visé (Bc proche de 1) ; +2. **une collection de paires spécialisées**, chacune sur un clade que les paires précédentes résolvent mal. + +### 20.3 Réponses de l'utilisateur + +1. **Granularité des clades à raffiner : elle dépend du nombre de paires autorisées dans le multiplex.** + Avec 2 paires : une paire générique + une seule paire spécialisée, placée assez haut dans la taxonomie, + autour d'un gros groupe très mal résolu (celle qui récupère le plus d'espèces). Avec 30 paires : on peut + éclater presque chaque famille ou chaque genre mal résolu. +2. **Éviter d'amplifier les autres clades : pas par défaut, mais une option.** Elle permet de concentrer + l'effort d'amplification sur le clade d'intérêt sans le disperser, surtout s'il est sous-représenté dans la + matrice d'ADN. +3. **Critère de résolution : une espèce est résolue si au moins un marqueur la résout.** +4. **Budget : deux formulations possibles, non tranchées** (maximiser les espèces résolues avec au plus K + paires ; ou minimiser le nombre de paires pour atteindre un taux donné). +5. **La paire de base est choisie d'abord**, ce qui élimine les espèces bien reconnues ; on raisonne ensuite + sur ce qui reste (glouton hiérarchique). + +### 20.4 Ce que cela implique (propositions, à valider) + +- **Algorithme glouton hiérarchique** : paire de base sur tout le groupe (tri Bc puis Bs, comme le mode A) ; + puis, tant que le budget le permet, parmi les nœuds de la taxonomie contenant des espèces non résolues, + choisir la paire spécialisée au **plus grand gain marginal** (espèces nouvellement résolues). La + granularité exigée en 20.3-1 en découle sans règle séparée : un petit budget favorise les nœuds hauts, un + grand budget finit par les genres. Les deux formulations du budget (20.3-4) sont le même gain marginal avec + deux conditions d'arrêt (K paires atteint, ou taux visé atteint) : le choix n'a pas à être tranché pour + écrire la boucle. +- **Réutilisation** : le mode A avec `-r ` produit déjà les paires conservées dans un clade et leur Bs. + Le mode A étant rapide (de 2 à 3 s pour Teleostei), faire une recherche par clade candidat est envisageable ; + le corpus devrait alors être chargé une seule fois et la sélection exemples/contre-exemples refaite en + mémoire pour chaque clade (aujourd'hui `LoadCorpus` relit les fichiers). +- **Manque 1 : les ensembles d'espèces résolues.** Aujourd'hui `PairStats` ne garde que des nombres + (`Identified`, `ExampleTaxa`). Il faut, pour chaque paire, l'ensemble des taxons identifiés et les classes + de taxons partageant un amplicon (même besoin qu'un portage de `-b`). +- **Question ouverte (moi) : relativement à quel ensemble une espèce est-elle « résolue » par un spécialiste + qui peut amplifier d'autres clades (20.3-2) ?** Le Bs actuel ne juge l'unicité de l'amplicon que parmi les + exemples. Si le spécialiste amplifie aussi des espèces extérieures au clade, un amplicon partagé avec un + extérieur ne permet pas d'identifier l'espèce dans un échantillon réel. Je propose de juger l'unicité sur + **tout ce que la paire amplifie dans la base de référence**. Cela demande un troisième rôle de séquence, + « contexte » : non exigée amplifiée, non pénalisée, mais comptée dans le partage d'amplicons. +- **Lien avec la voie B** : avec l'option « éviter les autres clades » (20.3-2), les voisins deviennent des + contre-exemples ; sans elle, ils sont du contexte. La voie B (amorces spécifiques d'une espèce ou de + quelques espèces proches) et l'étape de spécialisation ressemblent donc beaucoup : je propose de les + spécifier ensemble, comme un seul mécanisme avec une option d'exclusion. +- **Non abordé, à demander** : compatibilité des amorces d'un même multiplex (Tm voisines, absence de + dimères entre amorces de paires différentes). `obithermo` fournit déjà les Tm. +- **Sortie** : un ensemble ordonné de paires avec, pour chacune, les espèces qu'elle ajoute (gain marginal) ; + la courbe cumulée des espèces résolues en fonction du nombre de paires répond aux deux formulations du + budget. + +### 20.5 Prochaines étapes + +Pas de code avant accord sur 20.4. À trancher : (a) le jugement d'unicité (tout ce qui est amplifié, ou +seulement le clade), donc le rôle « contexte » ; (b) si les voies B et de spécialisation sont spécifiées +ensemble ; (c) les contraintes de multiplex. Le reste de l'inventaire du mode A (validation `-e 0`, nettoyage, +mémoire, tests shell, documentation) est indépendant et peut avancer. + +## 21. Suite de la spécification : résolution, voie B, optimisation de Tm (2026-11-03) + +**Statut : spécification, aucun code.** « Décision » = réponse de l'utilisateur ; « Proposition » et +« Question » = de moi, non validées. + +### 21.1 Décisions + +1. **Une espèce est résolue par rapport à tout ce qui est amplifié**, dans le clade cible **et** dans le clade + voisin non cible. Confirme le jugement d'unicité sur tout ce que la paire amplifie dans la référence. +2. **La voie B doit être spécifiée complètement** avant tout code (et avec la spécialisation, 20.4). +3. **Tm** : idéalement oui, des Tm voisines entre amorces. Aujourd'hui la longueur est fixée et la Tm subie. + Idée de l'utilisateur : **fixer la Tm et optimiser la longueur à la fin pour équilibrer les Tm**, et cela + **pour n'importe quelle paire**. Ajouter à `obiprimers` une **phase d'optimisation de Tm** avec deux + paramètres : **un écart de Tm maximal entre les deux amorces** et **une Tm cible pour les paires**. + +### 21.2 Conséquence de 21.1-1 : quatre rôles de séquences (proposition) + +| Rôle | Exigence sur l'amplification | Compte dans l'unicité des amplicons | +|---|---|---| +| Exemple | doit être amplifiée (quorum de sensibilité) | oui | +| Contre-exemple | ne doit pas l'être (quorum de faux positifs) | **oui, si amplifiée** | +| Contexte | libre, aucune pénalité | **oui, si amplifiée** | +| Ignorée | sans objet | non | + +Le rôle « contexte » est nouveau (aujourd'hui : exemple, contre-exemple, ignorée). L'option « éviter les autres +clades » (20.3-2) transforme simplement leurs séquences de contexte en contre-exemples. +**Écart avec ecoPrimers à noter** : le Bs d'ecoPrimers ne juge l'unicité que parmi les exemples, même quand des +contre-exemples sont amplifiés. Proposition : **garder le Bs d'ecoPrimers inchangé** (c'est ce qui reproduit +les tables de référence) et ajouter un second indice calculé sur tout ce qui est amplifié (nom à choisir, +p. ex. `bs_all`), qui sert à la sélection multi-marqueurs ; il coïncide avec `bs` quand il n'y a ni +contre-exemple amplifié ni contexte. + +### 21.3 Optimisation de Tm (proposition, paramètres de l'utilisateur) + +- **Paramètres** : `Tm cible` et `écart de Tm maximal`. Proposition d'interprétation : la cible porte sur la + **moyenne des deux amorces** et l'écart borne leur différence, ce qui revient à imposer à chaque amorce + `cible ± écart/2`. +- **Où** : phase **après** la sélection des paires, sur chaque paire retenue : choisir les longueurs + `(L1, L2)` des deux amorces satisfaisant la cible et l'écart. +- **Ce qui rend la phase possible** : raccourcir une amorce **par son extrémité 5'** ne peut que conserver ou + réduire ses erreurs sur chaque site (c'est une sous-chaîne), donc la sensibilité est conservée. Allonger + par le 5' demande de **revérifier** les sites. Raccourcir fait en revanche perdre de la spécificité + (l'amorce raccourcie peut s'apparier ailleurs), d'où la nécessité de **réévaluer la paire** (sensibilité, + faux positifs, Bc, Bs, Tm) après chaque changement de longueur. +- **Réutilisation** : `obithermo` (Tm, Tm minimale), `MatchWords` et `FindPairs` sur une poignée de mots : + le coût est celui d'une évaluation par variante de longueur et par paire. +- **Intérêt pour le multiplex** : une Tm cible commune à toutes les paires d'un ensemble les rend compatibles + par construction. +- **Questions** : (a) la cible porte-t-elle sur la Tm de l'appariement parfait ou sur la Tm minimale sur les + sites (la pire), ou sur une valeur intermédiaire ? (b) bornes de longueur [Lmin, Lmax] ? Le mot actuel est + limité à 31 bases (codage 2 bits sur 64). (c) extension et raccourcissement seulement par le 5', en gardant + l'ancrage 3' ? (d) « pour n'importe quelle paire » : phase intégrée à la recherche, ou aussi outil + indépendant appliqué à des amorces fournies par l'utilisateur ? (e) une paire dont la Tm ne peut pas être + satisfaite dans les bornes : écartée ou signalée ? + +### 21.4 Voie B : esquisse de spécification (proposition) + +**But.** Une paire qui amplifie les cibles T (une espèce, ou quelques espèces très proches) et, autant que +possible, aucun de leurs voisins non cibles C. Facteur de conception : l'espace de recherche est élagué en +se concentrant sur les régions **non conservées** (celles où T diffère de C), l'inverse du mode A. + +**Rôles.** T = exemples ; C = contre-exemples (voisins à ne pas amplifier) ; le reste, contexte ou ignoré. + +**Constat.** Le pipeline actuel fait déjà fonctionner un cas voie B de façon naïve (`-r 8028 -i 8015`, *Salmo* +contre les autres Salmonidae : 7 exemples, 76 contre-exemples, 15 534 mots stricts), mais sans élagage, et le +critère « un quorum de faux positifs de 10 % » y est trop permissif pour une amorce spécifique. + +**Principe d'une paire spécifique.** Une paire n'amplifie un voisin que si **ses deux amorces** s'y apparient. +Il suffit donc d'**une amorce spécifique** (loin de tous les sites de C) et d'une **amorce partenaire** +quelconque et conservée dans T, à distance d'amplicon permise. Les candidats se génèrent ainsi : mots de T +(quorum de T) dont aucun site à <= `e_exclusion` erreurs n'existe dans C (au-delà du quorum de faux +positifs) ; puis le partenaire dans la fenêtre d'amplicon de chaque occurrence dans T. + +**Deux tolérances (proposition).** `-e` (sensibilité sur T) et une tolérance d'exclusion **`e_exclusion`** sur C, +éventuellement plus grande : un voisin est « amplifié » s'il s'apparie à `e_exclusion` erreurs ou moins. Cela +donne une marge de sécurité sans dégrader la sensibilité sur T. + +**Score de marge (proposition).** Pour une amorce : nombre minimal d'erreurs d'un site de C ; pour la paire : +la meilleure des deux amorces (un voisin n'est amplifié que si les deux s'apparient), avec un poids +particulier pour les erreurs proches du 3'. Il compléterait Yule x Bs pour classer. + +**Sortie.** Colonnes du mode A plus la marge des deux amorces. + +**Échelle (à prévoir).** Pour T = un génome, le nombre de mots candidats est celui des positions du génome +(16 000 pour une mitochondrie, des millions pour une bactérie) : l'index de graines (`MatchWords`) indexe les +**mots** et balaie le texte de C ; pour des millions de mots il faudra les traiter par lots, chacun balayant C +une fois. À dimensionner selon l'usage visé. + +**Questions pour finaliser.** +1. **Exclusion stricte ou marge** : une paire est-elle admise seulement si elle n'amplifie aucun voisin à + `e_exclusion` erreurs, ou classe-t-on par marge sans seuil dur ? Faut-il une contrainte particulière sur + les 3' derniers nucléotides des voisins (au moins une erreur dans les 5 dernières bases) ? +2. **Tailles typiques** de T et de C (une espèce avec combien de séquences ? voisins = un genre, une famille, + toute la base de 8 000 génomes mitochondriaux, des génomes entiers de bactéries ?). +3. **Type de séquences** : génomes organellaires entiers (comme le tutoriel), loci courts, génomes complets ? +4. **Plusieurs cibles** (quelques espèces proches) : amplifier le groupe sans les distinguer (détection de + groupe), ou les amplifier et les **résoudre** entre elles (Bs au sein de T, comme la spécialisation) ? +5. **Une amorce spécifique + une générale** est-il acceptable, ou les deux amorces doivent-elles être + spécifiques ? +6. **Sensibilité sur T** : quorum comme en mode A (0,9), ou toutes les séquences de T (intra-spécifique) ? + +### 21.5 Travaux préalables communs (sans les engager) + +Rôle « contexte » ; indice de résolution sur tout ce qui est amplifié ; ensembles d'espèces résolues par +paire ; corpus chargé une seule fois et sélection refaite en mémoire ; réévaluation d'une paire pour des +longueurs d'amorces variables ; marge d'une amorce (erreurs minimales sur C). + +## 22. Réponses de l'utilisateur sur la résolution et la Tm (2026-11-03) + +**Statut : spécification, aucun code.** La voie B attend les réponses de l'utilisateur (« encore un peu »). + +### 22.1 Redéfinir Bs : oui (décision de principe, détails à confirmer) + +Question de l'utilisateur : ecoPrimers calcule Bs sur les seules séquences amplifiées parce que ce sont les +seules chargées en mémoire ; ne vaudrait-il pas mieux **charger tout le monde, chercher les amorces +seulement dans une sélection, et tester sur tout le monde** ? Réponse : **oui**, et cela redéfinit Bs. Trois +populations distinctes : + +| Population | Sert à | +|---|---| +| exemples (sélection de recherche) | générer les mots conservés ; sensibilité ; dénominateur de Bc et de Bs | +| contre-exemples | quorum de faux positifs | +| **tout le monde** (population d'évaluation) | juger l'unicité des amplicons : une espèce exemple est résolue si son amplicon n'est partagé avec aucun autre taxon amplifié, où qu'il soit dans la base | + +Proposition de définition : `Bs = espèces exemples dont l'amplicon est unique parmi tous les taxons amplifiés de +la population d'évaluation / espèces exemples amplifiées`. Elle coïncide avec le Bs d'ecoPrimers quand la +population d'évaluation est la sélection (aucun contre-exemple ni contexte amplifié). **Piège à ne pas +oublier** : avec `-i 1` (validation de la section 17) le Bs d'ecoPrimers ignore les amplicons des +contre-exemples amplifiés ; la nouvelle définition changera donc ces valeurs. Pour garder les tests de +non-régression, proposition : conserver aussi le Bs d'ecoPrimers dans une colonne à part (p. ex. +`bs_examples`) et définir `bs` comme le nouvel indice ; à confirmer. Colonne informative à ajouter : +nombre de taxons hors exemples amplifiés. + +**Mémoire (proposition)** : charger tout le monde n'est pas nécessaire pour un passage simple. Après la +recherche (sur la sélection, en mémoire), un **passage d'évaluation en flux** relit la base entière avec les +seules amorces des paires retenues (quelques centaines de mots : `MatchWords` en 0,04 s pour 65 Mb) et ne +conserve que, par paire, un hash d'amplicon et un taxon par séquence amplifiée, jamais les séquences. Coût : +une relecture du fichier (comme la double lecture de la taxonomie embarquée). Pour la sélection multi-marqueurs +(évaluations répétées par clade), une population d'évaluation **en mémoire sous forme compacte** (2 bits par +base, soit ~33 Mo pour les 131 Mb du tutoriel) serait préférable. À trancher à la spécification de la +sélection multi-marqueurs. + +### 22.2 Longueur des amorces : 61 nucléotides (question de l'utilisateur) + +L'utilisateur rappelle qu'agrep avait été réimplémenté en 64 bits pour des amorces jusqu'à 61 nucléotides +(`MAX_PAT_LEN` vaut 64 dans `obiapat/apat.h`). État **actuel** de obiprimers : **31 bases au maximum**, limite +inscrite dans `match.go:48`, `words.go:44`, `pairs.go:31`, dans `obikmer` (`encodekmer.go`, +`seqquorum.go:63` : k-mers sur 64 bits, 2 bits par base) et par des tampons fixes `[32]byte` +(`pairs_eval.go:208`). `obithermo` n'a pas de limite de longueur (tranches). + +La stratégie de l'index de graines **n'est pas bloquée** : elle se généralise à 61 bases (et 64) avec une fenêtre +de **128 bits** (deux mots de 64 bits, décalage et popcount sur les deux) ; les moitiés (<= 31 bases) tiennent +dans un mot et le chemin d'index haché pour les grandes graines est déjà utilisé et testé (k = 25 et 31). +L'oracle de test existant (`matchNaive`) fixe la justesse. Le compteur de mots stricts (k-mers sur 64 bits, +type partagé d'`obikmer`) resterait limité à 31 : **proposition** : chercher les mots conservés avec un k <= 31 +(le « coeur » conservé), puis obtenir les amorces plus longues par **extension** pendant la phase de Tm +(22.3), vérifiée par l'appariement généralisé. + +### 22.3 Optimisation de Tm : réponses aux questions de 21.3 + +(a) **Cible sur l'apparié non parfait** (la Tm de l'amorce sur les sites réels, avec leurs erreurs), pas sur +l'appariement parfait. *Question restante* : quelle statistique des Tm de sites ? Le minimum (`tm_min` d'ecoPrimers) +est dominé par des valeurs aberrantes : sur les tables de référence on lit des minimums de 12 à 45 °C pour des +amorces dont la Tm parfaite est de 52 à 62 °C ; une cible sur le minimum imposerait des longueurs absurdes. +Proposition : cible sur la **moyenne ou la médiane** des Tm de sites (éventuellement pondérées par taxon), le +minimum restant rapporté ; ou un quantile bas réglable. À choisir. +(b) bornes de longueur : traitées en 22.2 (jusqu'à 61 bases possible avec la généralisation). +(c) **Extension et raccourcissement du côté le mieux conservé** (et non seulement par le 5', contrairement à ma +proposition). Cela se calcule avec un profil de conservation par position (fraction des sites cibles qui +portent la même base, comme le PSSM du tutoriel) autour de chaque amorce, tiré des occurrences déjà +trouvées. Attention : la contrainte `-3` s'applique à l'extrémité 3' *de l'amorce*, quelle que soit sa +position sur la séquence ; étendre ou raccourcir du côté 3' déplace donc cette extrémité et modifie la +longueur d'amplicon. +(d) **Intégrée à la recherche, en phase de post-traitement** de la même exécution de obiprimers (pas d'outil +indépendant pour l'instant). +(e) Une paire dont la Tm n'est pas atteignable dans les bornes est **écartée**. + +### 22.4 En attente + +Réponses de l'utilisateur sur la voie B (questions de 21.4) ; choix de la statistique de Tm (22.3-a) ; nom et +sort de l'ancien Bs (22.1) ; en-têtes éventuels à prévoir pour les colonnes de population d'évaluation. + +## 23. Décisions : un seul Bs, tout en mémoire (2026-11-03) + +**Statut : spécification, aucun code.** Cette section **remplace** la proposition de 22.1 (colonne `bs_examples`, +passage d'évaluation en flux). + +### 23.1 Décisions de l'utilisateur + +1. **Un seul indice Bs, calculé d'une seule façon.** Pas de colonne `bs_examples`, pas de variante : « on décide + un indice, il se calcule d'une certaine façon, point barre ». Le Bs d'ecoPrimers disparaît. +2. **Tout charger en mémoire.** Ne pas tout charger interdirait de lire la base de référence sur l'entrée + standard, que les obitools permettent ; la base est lue par le mécanisme de lecture standard + (comme `obiconvert`), sans format dédié. Le passage d'évaluation en flux de 22.1 est abandonné. +3. **Définition de Bs** : le nombre de taxons **bien identifiés parmi les cibles**, en considérant **tous les + barcodes amplifiés, y compris ceux hors des cibles et ceux des contre-exemples**. Une espèce cible est + mal identifiée si son barcode, unique parmi les cibles, est partagé avec une espèce non cible ou un + contre-exemple. + +### 23.2 Définition formelle retenue (traduction de 23.1-3) + +Pour une paire, soit A l'ensemble des séquences qu'elle amplifie, **de n'importe quel rôle** (cibles, +contre-exemples, autres). Pour un amplicon a, T(a) est l'ensemble des taxons (au rang choisi) des séquences de A +qui le portent. Un taxon cible t est **bien identifié** s'il existe un amplicon a produit par t avec T(a) = {t}. + +### 23.3 Conséquences + +- **Fin de la notion de séquence « ignorée »** : les séquences qui ne sont ni cibles ni contre-exemples ne sont + plus écartées à la lecture (comportement actuel, hérité d'ecoPrimers) ; elles sont **non cibles** : elles ne + comptent ni dans les quorums de sensibilité et de faux positifs, mais leurs amplicons comptent dans Bs. + Pour exclure des séquences de l'analyse, on filtre l'entrée en amont (`obigrep`), à la manière des obitools. + Le rôle « contexte » de 21.2 devient donc simplement « non cible ». +- **Les valeurs changent, y compris pour le tutoriel** : une paire conservée chez tous les vertébrés amplifie + aussi des non-Cyprinidae ; son Bs baissera par rapport aux tables d'ecoPrimers, même sans `-i`. Les tests de + non-régression de la section 17 et de la section 18 ne peuvent plus comparer Bs/identifiés avec ecoPrimers sur + la base entière. **Proposition, sans option ni colonne supplémentaire** : au niveau de la bibliothèque, un corpus + construit avec les seules cibles et contre-exemples (c'est ce que fait `LoadCorpus` aujourd'hui) redonne + exactement la définition d'ecoPrimers, car il n'y a alors aucune séquence hors sélection à amplifier ; les tests de + non-régression continuent d'utiliser ce corpus, la CLI charge toujours tout. +- **Coût mémoire** : tout le fichier est conservé (131 Mb de bases pour 194 Mo de tas dans le tutoriel, contre + 5 Mb pour Cyprinidae aujourd'hui). Une représentation compacte (2 bits par base) reste possible plus tard ; + non décidée. +- **Coût de calcul** : l'appariement approché des mots et l'énumération des amplicons passent de la sélection à + tout le corpus (l'exécution `-r 7953 -i 1` du tutoriel, qui charge déjà tout, prend 1,9 s) ; les quorums ne + comptent toujours que cibles et contre-exemples. + +### 23.4 Points à confirmer (propositions, non validées) + +1. **Dénominateur de Bs** : taxons cibles **amplifiés** (comme le code d'ecoPrimers : 2 650 identifiés sur 3 805 + amplifiés = 0,696) ou **tous** les taxons cibles ? Avec le premier, Bc (couverture) et Bs (discrimination) + restent indépendants et `Bc x Bs` est la part des cibles résolues par la paire, ce qui est le gain mesuré par + la sélection multi-marqueurs. Avec le second, Bs <= Bc et ce produit disparaît. **Je recommande le + premier.** L'article (« nombre de taxons identifiés sur |E| ») et le code divergent déjà. +2. **Taxon à plusieurs amplicons** : bien identifié dès qu'**un** de ses amplicons est unique (règle d'ecoPrimers), + ou seulement si **tous** le sont ? Je propose de garder la règle d'ecoPrimers. +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) + +**Statut : spécification, aucun code.** Complète et corrige la section 23. + +### 24.1 Décisions de l'utilisateur + +1. **Dénominateur de Bs : les taxons cibles amplifiés** (comme ecoPrimers : 2 650 / 3 805 = 0,696). + Bc (couverture) et Bs (discrimination) restent indépendants ; `Bc x Bs` est la part des cibles résolues. +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 + 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.) + +### 24.2 ecoPrimers applique bien l'autre règle (vérifié) + +Lecture de `taxstats.c` : pour chaque amplicon distinct (parmi les exemples), on dresse l'ensemble de ses +taxons ; si cet ensemble n'a qu'**un** taxon (`taxoncount == 1`), ce taxon est marqué bien identifié. Un +taxon est donc compté **dès qu'un** de ses amplicons lui est propre. Confirmations : (i) mon implémentation de +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 = +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 +et les tables reproduites. + +### 24.3 Conséquences + +- **Les différences entre règles seront faibles sur le tutoriel** (un génome par espèce) et **plus fortes** + sur des bases à plusieurs séquences par espèce (variation intra-spécifique), typiquement des bases de loci courts. +- **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 + **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. +- **Implémentation (plus tard)** : la phase 2 de `FindPairs` doit enregistrer le hash d'amplicon de **toutes** les + séquences amplifiées (cibles, contre-exemples, non-cibles) et non plus des seuls exemples ; un taxon cible est + mal identifié si l'un de ses amplicons est porté par un autre taxon. + +### 24.4 Questions restantes + +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). diff --git a/pkg/obiprimer/search_real_test.go b/pkg/obiprimer/search_real_test.go index babb400..53b9466 100644 --- a/pkg/obiprimer/search_real_test.go +++ b/pkg/obiprimer/search_real_test.go @@ -16,6 +16,12 @@ import ( // loadTutorial reads the Cyprinidae examples of the tutorial data (the opt-in // tests need OBIPRIMERS_TUTORIAL, see TestStrictWordsAgainstEcoPrimers). func loadTutorial(t *testing.T) (*Corpus, *Taxa, *obitax.Taxonomy) { + t.Helper() + return loadTutorialWith(t, nil) +} + +// loadTutorialWith is loadTutorial with counterexample clades (ecoPrimers -i). +func loadTutorialWith(t *testing.T, ignore []string) (*Corpus, *Taxa, *obitax.Taxonomy) { t.Helper() dir := os.Getenv("OBIPRIMERS_TUTORIAL") if dir == "" { @@ -33,7 +39,11 @@ func loadTutorial(t *testing.T) (*Corpus, *Taxa, *obitax.Taxonomy) { if err != nil { t.Fatal(err) } - corpus, err := LoadCorpus(seqs, Selection{Restrict: restrict}) + sel := Selection{Restrict: restrict} + if sel.Ignore, err = BelongingTo(taxonomy, ignore); err != nil { + t.Fatal(err) + } + corpus, err := LoadCorpus(seqs, sel) if err != nil { t.Fatal(err) } @@ -173,3 +183,61 @@ func TestReferenceSitesAreRealSites(t *testing.T) { }) } } + +// ecoPrimers with counterexamples: -r 7953 -i 1 (the Cyprinidae against all the +// other sequences of the database, 7497 counterexamples of 7387 species), -e 3 +// -3 2 -l 30 -L 150 -c. ecoPrimers needed 45 minutes and keeps 173 pairs. +// testdata/ecoprimers_Cyprinidae_I1.pairs.tsv holds for each: the primers, the +// good flags, the amplified example and counterexample sequences, the amplified +// example and counterexample taxa, and the identified taxa. +func TestCounterExamplesAgainstEcoPrimers(t *testing.T) { + corpus, taxa, _ := loadTutorialWith(t, []string{"1"}) + 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)) + if err != nil { + t.Fatal(err) + } + t.Logf("%+v", *info) + + f, err := os.Open(filepath.Join("testdata", "ecoprimers_Cyprinidae_I1.pairs.tsv")) + if err != nil { + t.Fatal(err) + } + defer f.Close() + type row struct { + good string + v [5]int + } + want := map[string]row{} + sc := bufio.NewScanner(f) + for sc.Scan() { + var p1, p2 string + var r row + if _, err := fmt.Sscanf(sc.Text(), "%s\t%s\t%s\t%d\t%d\t%d\t%d\t%d", &p1, &p2, &r.good, + &r.v[0], &r.v[1], &r.v[2], &r.v[3], &r.v[4]); err != nil { + t.Fatal(err) + } + want[p1+"/"+p2] = r + } + + if len(pairs) != len(want) { + t.Errorf("%d pairs, ecoPrimers reports %d", len(pairs), len(want)) + } + bad := 0 + for _, g := range pairs { + w, ok := want[g.Primer1+"/"+g.Primer2] + if !ok { + t.Errorf("pair %s/%s is not in the ecoPrimers result", g.Primer1, g.Primer2) + 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 { + 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) + } + } + } +} diff --git a/pkg/obiprimer/testdata/ecoprimers_Cyprinidae_I1.pairs.tsv b/pkg/obiprimer/testdata/ecoprimers_Cyprinidae_I1.pairs.tsv new file mode 100644 index 0000000..ae1da51 --- /dev/null +++ b/pkg/obiprimer/testdata/ecoprimers_Cyprinidae_I1.pairs.tsv @@ -0,0 +1,173 @@ +AGCAATCCAGGTCAGTTT TGCCACCCCAGTATGTCC bG 309 454 301 452 197 +GCAATCCAGGTCAGTTTC TGCCACCCCAGTATGTCC bG 309 455 301 453 197 +AGCAATCCAGGTCAGTTT CCACCCCAGTATGTCCTT bG 311 484 302 482 197 +CCACCCCAGTATGTCCTT GCAATCCAGGTCAGTTTC Gb 311 485 302 483 197 +AGCAATCCAGGTCAGTTT GCCACCCCAGTATGTCCT bG 308 445 300 444 196 +GCAATCCAGGTCAGTTTC GCCACCCCAGTATGTCCT bG 308 446 300 445 196 +ACCGGAGCAATCCAGGTC TGCCACCCCAGTATGTCC bG 308 453 300 451 196 +GAGCAATCCAGGTCAGTT TGCCACCCCAGTATGTCC bG 308 453 300 451 196 +AGACCGGAGCAATCCAGG CCACCCCAGTATGTCCTT bG 310 483 301 481 196 +CCACCCCAGTATGTCCTT GAGCAATCCAGGTCAGTT Gb 310 483 301 481 196 +ACCGGAGCAATCCAGGTC CCACCCCAGTATGTCCTT bG 310 484 301 482 196 +CCACCCCAGTATGTCCTT GACCGGAGCAATCCAGGT Gb 310 484 301 482 196 +CCACCCCAGTATGTCCTT GGTCAGTTTCTATCTGTA Gb 309 460 300 459 195 +AGGTCAGTTTCTATCTGT TGCCACCCCAGTATGTCC bG 307 438 299 437 195 +CAGGTCAGTTTCTATCTG TGCCACCCCAGTATGTCC bG 307 440 299 439 195 +GGTCAGTTTCTATCTGTA TGCCACCCCAGTATGTCC bG 307 441 299 440 195 +GAGCAATCCAGGTCAGTT GCCACCCCAGTATGTCCT bG 307 444 299 443 195 +ACCGGAGCAATCCAGGTC GCCACCCCAGTATGTCCT bG 307 445 299 444 195 +GACCGGAGCAATCCAGGT GCCACCCCAGTATGTCCT bG 307 445 299 444 195 +AGGTCAGTTTCTATCTGT CCACCCCAGTATGTCCTT bG 309 469 300 467 195 +CAGGTCAGTTTCTATCTG CCACCCCAGTATGTCCTT bG 309 470 300 468 195 +GGAGCAATCCAGGTCAGT TGCCACCCCAGTATGTCC bG 307 454 299 452 195 +CCAGGTCAGTTTCTATCT TGCCACCCCAGTATGTCC bG 307 455 299 453 195 +AATCCAGGTCAGTTTCTA TGCCACCCCAGTATGTCC bG 307 455 299 453 195 +CAATCCAGGTCAGTTTCT TGCCACCCCAGTATGTCC bG 307 455 299 453 195 +ATCCAGGTCAGTTTCTAT TGCCACCCCAGTATGTCC bG 307 455 299 453 195 +TCCAGGTCAGTTTCTATC TGCCACCCCAGTATGTCC bG 307 455 299 453 195 +CCACCCCAGTATGTCCTT GGAGCAATCCAGGTCAGT Gb 309 484 300 482 195 +CAATCCAGGTCAGTTTCT CCACCCCAGTATGTCCTT bG 309 485 300 483 195 +AATCCAGGTCAGTTTCTA CCACCCCAGTATGTCCTT bG 309 485 300 483 195 +ATCCAGGTCAGTTTCTAT CCACCCCAGTATGTCCTT bG 309 485 300 483 195 +CCACCCCAGTATGTCCTT TCCAGGTCAGTTTCTATC Gb 309 485 300 483 195 +CCACCCCAGTATGTCCTT CCAGGTCAGTTTCTATCT Gb 309 485 300 483 195 +GCCACCCCAGTATGTCCT GGTCAGTTTCTATCTGTA Gb 306 437 298 437 194 +AGGTCAGTTTCTATCTGT GCCACCCCAGTATGTCCT bG 306 440 298 439 194 +CAGGTCAGTTTCTATCTG GCCACCCCAGTATGTCCT bG 306 442 298 441 194 +GCCACCCCAGTATGTCCT GGAGCAATCCAGGTCAGT Gb 306 445 298 444 194 +GCCACCCCAGTATGTCCT TCCAGGTCAGTTTCTATC Gb 306 446 298 445 194 +CCAGGTCAGTTTCTATCT GCCACCCCAGTATGTCCT bG 306 446 298 445 194 +CAATCCAGGTCAGTTTCT GCCACCCCAGTATGTCCT bG 306 446 298 445 194 +AATCCAGGTCAGTTTCTA GCCACCCCAGTATGTCCT bG 306 446 298 445 194 +ATCCAGGTCAGTTTCTAT GCCACCCCAGTATGTCCT bG 306 446 298 445 194 +CGGAGCAATCCAGGTCAG TGCCACCCCAGTATGTCC bG 304 445 296 444 192 +CCGGAGCAATCCAGGTCA TGCCACCCCAGTATGTCC bG 304 445 296 444 192 +TGCCACCCCAGTATGTCC TTTTCCTAGTACGAAAGG Gb 308 453 300 452 192 +CCACCCCAGTATGTCCTT CCGGAGCAATCCAGGTCA Gb 306 476 297 474 192 +CCACCCCAGTATGTCCTT CGGAGCAATCCAGGTCAG Gb 306 476 297 474 192 +GACCGGAGCAATCCAGGT TGCCACCCCAGTATGTCC bG 304 453 296 451 192 +CCACCCCAGTATGTCCTT GTCAGTTTCTATCTGTAA Gb 304 408 295 407 190 +CCACCCCAGTATGTCCTT TTTTCCTAGTACGAAAGG Gb 310 484 301 482 192 +CAGACCGGAGCAATCCAG CCACCCCAGTATGTCCTT bG 306 484 297 482 192 +GTCAGTTTCTATCTGTAA TGCCACCCCAGTATGTCC bG 302 390 294 390 190 +AGCAATCCAGGTCAGTTT CCCAGTATGTCCTTATTT bG 301 310 292 309 187 +CCCAGTATGTCCTTATTT GCAATCCAGGTCAGTTTC Gb 301 311 292 310 187 +ACTTTTCCTAGTACGAAA TGCCACCCCAGTATGTCC bG 307 443 299 442 191 +GCCACCCCAGTATGTCCT TTTTCCTAGTACGAAAGG Gb 307 445 299 444 191 +CGGAGCAATCCAGGTCAG GCCACCCCAGTATGTCCT bG 303 443 295 442 191 +CCGGAGCAATCCAGGTCA GCCACCCCAGTATGTCCT bG 303 443 295 442 191 +CTTTTCCTAGTACGAAAG TGCCACCCCAGTATGTCC bG 307 448 299 447 191 +AGACCGGAGCAATCCAGG GCCACCCCAGTATGTCCT bG 303 445 295 444 191 +ACTTTTCCTAGTACGAAA CCACCCCAGTATGTCCTT bG 309 478 300 476 191 +CCACCCCAGTATGTCCTT CTTTTCCTAGTACGAAAG Gb 309 482 300 480 191 +GCCACCCCAGTATGTCCT GTCAGTTTCTATCTGTAA Gb 301 390 293 390 189 +AGACCGGAGCAATCCAGG CCCAGTATGTCCTTATTT bG 300 310 291 309 186 +CAGACCGGAGCAATCCAG CCCAGTATGTCCTTATTT bG 300 310 291 309 186 +CCCAGTATGTCCTTATTT GAGCAATCCAGGTCAGTT Gb 300 310 291 309 186 +ACCGGAGCAATCCAGGTC CCCAGTATGTCCTTATTT bG 300 311 291 310 186 +CCCAGTATGTCCTTATTT GACCGGAGCAATCCAGGT Gb 300 311 291 310 186 +CCCAGTATGTCCTTATTT TCAGACCGGAGCAATCCA Gb 300 311 291 310 186 +CCCAGTATGTCCTTATTT TTCAGACCGGAGCAATCC Gb 300 311 291 310 186 +CCCAGTATGTCCTTATTT GTTCAGACCGGAGCAATC Gb 300 311 291 310 186 +CTTTTCCTAGTACGAAAG GCCACCCCAGTATGTCCT bG 306 443 298 442 190 +ACTTTTCCTAGTACGAAA GCCACCCCAGTATGTCCT bG 306 443 298 442 190 +AGACCGGAGCAATCCAGG TGCCACCCCAGTATGTCC bG 302 453 294 451 190 +CCACCCCAGTATGTCCTT TCAGACCGGAGCAATCCA Gb 304 485 295 483 190 +AGGTCAGTTTCTATCTGT CCCAGTATGTCCTTATTT bG 299 309 290 308 185 +CAGGTCAGTTTCTATCTG CCCAGTATGTCCTTATTT bG 299 309 290 308 185 +CCCAGTATGTCCTTATTT GGTCAGTTTCTATCTGTA Gb 299 309 290 308 185 +ATCCAGGTCAGTTTCTAT CCCAGTATGTCCTTATTT bG 299 311 290 310 185 +CCAGGTCAGTTTCTATCT CCCAGTATGTCCTTATTT bG 299 311 290 310 185 +CCCAGTATGTCCTTATTT GGAGCAATCCAGGTCAGT Gb 299 311 290 310 185 +CCCAGTATGTCCTTATTT TCCAGGTCAGTTTCTATC Gb 299 311 290 310 185 +CAATCCAGGTCAGTTTCT CCCAGTATGTCCTTATTT bG 299 311 290 310 185 +AATCCAGGTCAGTTTCTA CCCAGTATGTCCTTATTT bG 299 311 290 310 185 +TGCCACCCCAGTATGTCC TTTCCTAGTACGAAAGGA Gb 305 449 297 447 189 +CAGACCGGAGCAATCCAG GCCACCCCAGTATGTCCT bG 301 446 293 445 189 +CCACCCCAGTATGTCCTT TTTCCTAGTACGAAAGGA Gb 307 479 298 477 189 +CCACCCCAGTATGTCCTT TACGAAAGGATCGGAAAA Gb 307 453 298 452 188 +CCCAGTATGTCCTTATTT CCGGAGCAATCCAGGTCA Gb 298 309 289 308 184 +CCCAGTATGTCCTTATTT CGGAGCAATCCAGGTCAG Gb 298 309 289 308 184 +GCCACCCCAGTATGTCCT TTTCCTAGTACGAAAGGA Gb 304 440 296 439 188 +GCAATTTACCATGCTCTG TTTTCCTAGTACGAAAGG Gb 308 505 299 501 189 +ACTTTTCCTAGTACGAAA GCAATTTACCATGCTCTG bG 307 478 298 474 188 +CAGACCGGAGCAATCCAG TGCCACCCCAGTATGTCC bG 300 454 292 452 188 +CCTAGTACGAAAGGATCG TGCCACCCCAGTATGTCC bG 303 427 295 425 187 +CCACCCCAGTATGTCCTT TTCAGACCGGAGCAATCC Gb 302 485 293 483 188 +TACGAAAGGATCGGAAAA TGCCACCCCAGTATGTCC bG 304 430 296 430 187 +GCCACCCCAGTATGTCCT TACGAAAGGATCGGAAAA Gb 304 430 296 430 187 +CCACCCCAGTATGTCCTT CCTAGTACGAAAGGATCG Gb 305 459 296 457 187 +CTTTTCCTAGTACGAAAG GCAATTTACCATGCTCTG bG 307 498 298 494 188 +CGCAATTTACCATGCTCT TTTTCCTAGTACGAAAGG Gb 306 463 297 459 187 +GCAATTTACCATGCTCTG TTTCCTAGTACGAAAGGA Gb 307 506 298 502 188 +GCAATTTACCATGCTCTG GGTCAGTTTCTATCTGTA Gb 301 433 292 429 186 +GCCACCCCAGTATGTCCT TCAGACCGGAGCAATCCA Gb 299 446 291 445 187 +CCTAGTACGAAAGGATCG GCCACCCCAGTATGTCCT bG 302 419 294 418 186 +CCCAGTATGTCCTTATTT GTCAGTTTCTATCTGTAA Gb 295 264 286 263 181 +ACTTTTCCTAGTACGAAA CGCAATTTACCATGCTCT bG 305 449 296 445 186 +TCGCAATTTACCATGCTC TTTTCCTAGTACGAAAGG Gb 305 456 296 452 186 +TTCGCAATTTACCATGCT TTTTCCTAGTACGAAAGG Gb 305 456 296 452 186 +AGTTCAGACCGGAGCAAT CCCAGTATGTCCTTATTT bG 296 310 287 309 182 +CGCAATTTACCATGCTCT CTTTTCCTAGTACGAAAG Gb 305 460 296 456 186 +CGCAATTTACCATGCTCT TTTCCTAGTACGAAAGGA Gb 305 463 296 459 186 +CCTAGTACGAAAGGATCG GCAATTTACCATGCTCTG bG 305 472 296 468 186 +TCAGACCGGAGCAATCCA TGCCACCCCAGTATGTCC bG 298 454 290 452 186 +ACTTTTCCTAGTACGAAA TCGCAATTTACCATGCTC bG 304 447 295 443 185 +ACTTTTCCTAGTACGAAA TTCGCAATTTACCATGCT bG 304 447 295 443 185 +CTTTTCCTAGTACGAAAG TCGCAATTTACCATGCTC bG 304 453 295 449 185 +CTTTTCCTAGTACGAAAG TTCGCAATTTACCATGCT bG 304 453 295 449 185 +CCACCCCAGTATGTCCTT GTTCAGACCGGAGCAATC Gb 300 485 291 483 186 +TCGCAATTTACCATGCTC TTTCCTAGTACGAAAGGA Gb 304 456 295 452 185 +TTTCGCAATTTACCATGC TTTTCCTAGTACGAAAGG Gb 304 456 295 452 185 +TTCGCAATTTACCATGCT TTTCCTAGTACGAAAGGA Gb 304 456 295 452 185 +GCAATTTACCATGCTCTG TACGAAAGGATCGGAAAA Gb 304 427 295 425 184 +CCTAGTACGAAAGGATCG CGCAATTTACCATGCTCT bG 303 433 294 429 184 +AGGTCAGTTTCTATCTGT GCAATTTACCATGCTCTG bG 299 431 290 427 184 +GCCACCCCAGTATGTCCT TTCAGACCGGAGCAATCC Gb 297 446 289 445 185 +ACTTTTCCTAGTACGAAA TTTCGCAATTTACCATGC bG 303 447 294 443 184 +GCAATTTACCATGCTCTG GTCAGTTTCTATCTGTAA Gb 295 375 286 373 182 +CTTTTCCTAGTACGAAAG TTTCGCAATTTACCATGC bG 303 453 294 449 184 +TTTCCTAGTACGAAAGGA TTTCGCAATTTACCATGC bG 303 456 294 452 184 +CCCAGTATGTCCTTATTT TTTTCCTAGTACGAAAGG Gb 300 311 291 310 180 +CCCAGTATGTCCTTATTT GAGTTCAGACCGGAGCAA Gb 294 310 285 309 180 +CCTAGTACGAAAGGATCG TCGCAATTTACCATGCTC bG 302 430 293 426 183 +CCTAGTACGAAAGGATCG TTCGCAATTTACCATGCT bG 302 430 293 426 183 +CGCAATTTACCATGCTCT TACGAAAGGATCGGAAAA Gb 302 420 293 418 182 +ACTTTTCCTAGTACGAAA CCCAGTATGTCCTTATTT bG 299 311 290 310 179 +CCCAGTATGTCCTTATTT CTTTTCCTAGTACGAAAG Gb 299 311 290 310 179 +CGCAATTTACCATGCTCT GGTCAGTTTCTATCTGTA Gb 297 426 288 422 182 +CCTAGTACGAAAGGATCG TTTCGCAATTTACCATGC bG 301 430 292 426 182 +CAGGTCAGTTTCTATCTG GCAATTTACCATGCTCTG bG 297 431 288 427 182 +TACGAAAGGATCGGAAAA TTCGCAATTTACCATGCT bG 301 418 292 416 181 +TACGAAAGGATCGGAAAA TCGCAATTTACCATGCTC bG 301 418 292 416 181 +GGCCTTTCGCAATTTACC TTTTCCTAGTACGAAAGG Gb 302 510 293 506 183 +ACTTTTCCTAGTACGAAA GGCCTTTCGCAATTTACC bG 301 479 292 475 182 +TACGAAAGGATCGGAAAA TTTCGCAATTTACCATGC bG 300 418 291 416 180 +CCCAGTATGTCCTTATTT TTTCCTAGTACGAAAGGA Gb 297 310 288 309 177 +CTTTTCCTAGTACGAAAG GGCCTTTCGCAATTTACC bG 301 500 292 496 182 +AGGTCAGTTTCTATCTGT CGCAATTTACCATGCTCT bG 295 428 286 424 180 +GGCCTTTCGCAATTTACC TTTCCTAGTACGAAAGGA Gb 301 511 292 507 182 +GCCTTTCGCAATTTACCA TTTTCCTAGTACGAAAGG Gb 299 440 290 436 180 +CTTTCGCAATTTACCATG TTTTCCTAGTACGAAAGG Gb 298 412 289 408 179 +CCCAGTATGTCCTTATTT CCTAGTACGAAAGGATCG Gb 296 305 287 304 176 +CCCAGTATGTCCTTATTT TACGAAAGGATCGGAAAA Gb 297 307 288 306 176 +GGTCAGTTTCTATCTGTA TCGCAATTTACCATGCTC bG 294 425 285 421 179 +ACTTTTCCTAGTACGAAA GCCTTTCGCAATTTACCA bG 298 436 289 432 179 +CTTTTCCTAGTACGAAAG GCCTTTCGCAATTTACCA bG 298 438 289 434 179 +GCCTTTCGCAATTTACCA TTTCCTAGTACGAAAGGA Gb 298 440 289 436 179 +CCTAGTACGAAAGGATCG GGCCTTTCGCAATTTACC bG 299 484 290 480 180 +ACTTTTCCTAGTACGAAA CTTTCGCAATTTACCATG bG 297 411 288 407 178 +CTTTCGCAATTTACCATG CTTTTCCTAGTACGAAAG Gb 297 411 288 407 178 +CTTTCGCAATTTACCATG TTTCCTAGTACGAAAGGA Gb 297 412 288 408 178 +GGCCTTTCGCAATTTACC TACGAAAGGATCGGAAAA Gb 298 425 289 423 178 +AGTTCAGACCGGAGCAAT CCACCCCAGTATGTCCTT bG 293 484 284 482 179 +CCTAGTACGAAAGGATCG GCCTTTCGCAATTTACCA bG 296 416 287 412 177 +CCTAGTACGAAAGGATCG CTTTCGCAATTTACCATG bG 295 390 286 386 176 +CCAGGTCAGTTTCTATCT GCAATTTACCATGCTCTG bG 294 513 285 509 179 +CCTTTCGCAATTTACCAT TTTTCCTAGTACGAAAGG Gb 295 411 286 407 176 +GCAATTTACCATGCTCTG TTCCTAGTACGAAAGGAT Gb 293 423 284 419 176 +CTTTCGCAATTTACCATG TACGAAAGGATCGGAAAA Gb 295 398 286 396 175 +ACTTTTCCTAGTACGAAA CCTTTCGCAATTTACCAT bG 294 409 285 405 175 +CCTTTCGCAATTTACCAT CTTTTCCTAGTACGAAAG Gb 294 409 285 405 175 +CCTTTCGCAATTTACCAT TTTCCTAGTACGAAAGGA Gb 294 410 285 406 175 +GCCTTTCGCAATTTACCA TACGAAAGGATCGGAAAA Gb 295 413 286 411 175 diff --git a/pkg/obiseq/taxonomy_methods.go b/pkg/obiseq/taxonomy_methods.go index e11356d..4934bcb 100644 --- a/pkg/obiseq/taxonomy_methods.go +++ b/pkg/obiseq/taxonomy_methods.go @@ -176,6 +176,10 @@ func (sequence *BioSequence) SetFamily(taxonomy *obitax.Taxonomy) *obitax.Taxon return sequence.SetTaxonAtRank(taxonomy, "family") } +// TaxonomicPathAttribute is the attribute holding the taxonomic path of a +// sequence, from the root to its taxon (see obiannotate --taxonomic-path). +const TaxonomicPathAttribute = "taxonomic_path" + func (sequence *BioSequence) SetPath(taxonomy *obitax.Taxonomy) []string { taxon := sequence.Taxon(taxonomy) path := taxon.Path() @@ -186,13 +190,13 @@ func (sequence *BioSequence) SetPath(taxonomy *obitax.Taxonomy) []string { spath[lpath-i] = path.Get(i).FullString(taxonomy.Code()) } - sequence.SetAttribute("taxonomic_path", spath) + sequence.SetAttribute(TaxonomicPathAttribute, spath) return spath } func (sequence *BioSequence) Path() []string { - path, ok := sequence.GetAttribute("taxonomic_path") + path, ok := sequence.GetAttribute(TaxonomicPathAttribute) if !ok { if taxo := obitax.DefaultTaxonomy(); taxo != nil { diff --git a/pkg/obitools/obiprimers/embedded_test.go b/pkg/obitools/obiprimers/embedded_test.go index daf0b3e..8f6f7ed 100644 --- a/pkg/obitools/obiprimers/embedded_test.go +++ b/pkg/obitools/obiprimers/embedded_test.go @@ -139,3 +139,59 @@ func TestPrimerPairsNeedsATaxonomy(t *testing.T) { t.Error("a nil taxonomy must be refused") } } + +// hasPath reports whether any sequence still carries its taxonomic path, and +// how many still have their taxid (which the search needs). +func scanSequences(it obiiter.IBioSequence) (withPath, withTaxid, total int) { + for it.Next() { + batch := it.Get() + for _, s := range batch.Slice() { + total++ + if s.HasAttribute(obiseq.TaxonomicPathAttribute) { + withPath++ + } + if s.HasAttribute("taxid") { + withTaxid++ + } + } + } + return +} + +// Once the taxonomy is built the paths only take memory: they are dropped, +// whatever the way the taxonomy was obtained, and nothing else is. +func TestTaxonomicPathIsDroppedAfterExtraction(t *testing.T) { + path := writeTemp(t, embeddedFasta) + + sequences, _, err := readSequencesAndTaxonomy([]string{path}, nil) // two passes over a file + if err != nil { + t.Fatal(err) + } + if withPath, withTaxid, total := scanSequences(sequences); total != 3 || withPath != 0 || withTaxid != 3 { + t.Errorf("file: %d sequences, %d with a path left, %d with a taxid", total, withPath, withTaxid) + } + + other := obitax.NewTaxonomy("other", "other", obitax.DefaultTaxidAlphabet) // -t given + sequences, _, err = readSequencesAndTaxonomy([]string{path}, other) + if err != nil { + t.Fatal(err) + } + if withPath, _, total := scanSequences(sequences); total != 3 || withPath != 0 { + t.Errorf("with -t: %d sequences, %d with a path left", total, withPath) + } + + var slice obiseq.BioSequenceSlice // standard input + for _, id := range []string{"a", "b"} { + seq := obiseq.NewBioSequence(id, []byte("acgtacgt"), "") + seq.SetAttribute("taxid", "taxon:10 [X]@species") + seq.SetAttribute(obiseq.TaxonomicPathAttribute, []string{"taxon:1 [root]@no rank", "taxon:10 [X]@species"}) + slice = append(slice, seq) + } + sequences, _, err = embeddedTaxonomyInMemory(obiiter.IBatchOver("stdin", slice, 10)) + if err != nil { + t.Fatal(err) + } + if withPath, withTaxid, total := scanSequences(sequences); total != 2 || withPath != 0 || withTaxid != 2 { + t.Errorf("stdin: %d sequences, %d with a path left, %d with a taxid", total, withPath, withTaxid) + } +} diff --git a/pkg/obitools/obiprimers/obiprimers.go b/pkg/obitools/obiprimers/obiprimers.go index 84f5a7a..21d036e 100644 --- a/pkg/obitools/obiprimers/obiprimers.go +++ b/pkg/obitools/obiprimers/obiprimers.go @@ -6,6 +6,7 @@ import ( "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault" "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiiter" "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiprimer" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitax" "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obiconvert" log "github.com/sirupsen/logrus" @@ -14,6 +15,16 @@ import ( const noTaxonomyAdvice = "no taxonomy: use the -t option, or give sequences annotated with their " + "taxonomic path (obiannotate --taxonomic-path)" +// dropTaxonomicPath deletes the taxonomic path of the sequences. The primer +// search needs the taxa of the sequences, which it reads from their taxid, not +// their path: once the taxonomy is built, the paths (a list of about twenty +// strings per sequence) only take memory. +func dropTaxonomicPath(sequences obiiter.IBioSequence) obiiter.IBioSequence { + return sequences.MakeIWorker(obiseq.AnnotatorToSeqWorker(func(s *obiseq.BioSequence) { + s.DeleteAttribute(obiseq.TaxonomicPathAttribute) + }), false) +} + // embeddedTaxonomyInMemory extracts the taxonomy from the taxonomic paths // carried by the sequences, which it has to hold in memory since the iterator // cannot be read twice. It returns an iterator over the same sequences. @@ -23,6 +34,10 @@ func embeddedTaxonomyInMemory(sequences obiiter.IBioSequence) (obiiter.IBioSeque if err != nil { return obiiter.NilIBioSequence, nil, fmt.Errorf("%s: %w", noTaxonomyAdvice, err) } + // every sequence is held in memory here: free the paths right away + for _, s := range slice { + s.DeleteAttribute(obiseq.TaxonomicPathAttribute) + } return obiiter.IBatchOver(source, slice, obidefault.BatchSize()), taxonomy, nil } @@ -39,7 +54,10 @@ func embeddedTaxonomyInMemory(sequences obiiter.IBioSequence) (obiiter.IBioSeque func readSequencesAndTaxonomy(args []string, given *obitax.Taxonomy) (obiiter.IBioSequence, *obitax.Taxonomy, error) { if given != nil { sequences, err := obiconvert.CLIReadBioSequences(args...) - return sequences, given, err + if err != nil { + return obiiter.NilIBioSequence, nil, err + } + return dropTaxonomicPath(sequences), given, nil } if len(args) == 0 { @@ -66,7 +84,10 @@ func readSequencesAndTaxonomy(args []string, given *obitax.Taxonomy) (obiiter.IB log.Infof("Embedded taxonomy extracted: %d taxa", taxonomy.Len()) sequences, err := obiconvert.CLIReadBioSequences(args...) - return sequences, taxonomy, err + if err != nil { + return obiiter.NilIBioSequence, nil, err + } + return dropTaxonomicPath(sequences), taxonomy, nil } // CLIReadSequencesAndTaxonomy reads the sequences of the command line and diff --git a/pkg/obiutils/set.go b/pkg/obiutils/set.go index c6c71b3..d73ad90 100644 --- a/pkg/obiutils/set.go +++ b/pkg/obiutils/set.go @@ -1,6 +1,10 @@ package obiutils -import "fmt" +import ( + "fmt" + "reflect" + "slices" +) type Set[E comparable] map[E]struct{} @@ -61,12 +65,43 @@ func (s Set[E]) Members() []E { return result } +// compareMembers orders two set members: numbers by value, strings +// alphabetically, anything else by its printed form. It is only used to give +// the string representation of a set a stable order. +func compareMembers(a, b any) int { + va, vb := reflect.ValueOf(a), reflect.ValueOf(b) + switch { + case va.CanInt() && vb.CanInt(): + return cmpOrdered(va.Int(), vb.Int()) + case va.CanUint() && vb.CanUint(): + return cmpOrdered(va.Uint(), vb.Uint()) + case va.CanFloat() && vb.CanFloat(): + return cmpOrdered(va.Float(), vb.Float()) + case va.Kind() == reflect.String && vb.Kind() == reflect.String: + return cmpOrdered(va.String(), vb.String()) + } + return cmpOrdered(fmt.Sprint(a), fmt.Sprint(b)) +} + +func cmpOrdered[T int64 | uint64 | float64 | string](a, b T) int { + switch { + case a < b: + return -1 + case a > b: + return 1 + } + return 0 +} + // String returns a string representation of the set. // -// It returns a string representation of the set by formatting the set's members using the fmt.Sprintf function. -// The resulting string is then returned. +// The members are formatted with fmt, in increasing order: a Go map has no +// order, and a representation changing from one call to the next is of no use +// to print, log or compare. func (s Set[E]) String() string { - return fmt.Sprintf("%v", s.Members()) + members := s.Members() + slices.SortFunc(members, func(a, b E) int { return compareMembers(a, b) }) + return fmt.Sprintf("%v", members) } // Union returns a new set that is the union of the current set and the specified set. diff --git a/pkg/obiutils/set_test.go b/pkg/obiutils/set_test.go index 1754bbb..fa6688c 100644 --- a/pkg/obiutils/set_test.go +++ b/pkg/obiutils/set_test.go @@ -300,3 +300,29 @@ func TestIntersection(t *testing.T) { t.Errorf("Intersection of two sets with no common elements returned incorrect result. Expected: %v, got: %v", expectedResult4, result4) } } + +// The representation of a set must not depend on the (random) order of the map. +func TestSetStringIsStable(t *testing.T) { + ints := NewSet(30, 4, 100, 2, 7, 15, 1) + for i := 0; i < 50; i++ { + if got := ints.String(); got != "[1 2 4 7 15 30 100]" { + t.Fatalf("ints: %s", got) + } + } + if got := NewSet("pear", "apple", "fig").String(); got != "[apple fig pear]" { + t.Errorf("strings: %s", got) + } + if got := NewSet(2.5, -1.0, 0.25).String(); got != "[-1 0.25 2.5]" { + t.Errorf("floats: %s", got) + } + type id int + if got := NewSet(id(3), id(1), id(2)).String(); got != "[1 2 3]" { + t.Errorf("named integer type: %s", got) + } + type pair struct{ a, b int } + for i := 0; i < 20; i++ { // no natural order: still stable + if got := NewSet(pair{2, 1}, pair{1, 2}, pair{1, 1}).String(); got != "[{1 1} {1 2} {2 1}]" { + t.Fatalf("structs: %s", got) + } + } +}