diff --git a/OBIPRIMERS_NOTES.md b/OBIPRIMERS_NOTES.md index 76fc038..20b1712 100644 --- a/OBIPRIMERS_NOTES.md +++ b/OBIPRIMERS_NOTES.md @@ -410,3 +410,66 @@ paires 1,1 s) sont celles à profiler si l'on veut aller plus loin ; le tri des dans `SeqKmerCounter` et la phase 2 de `FindPairs` en sont les candidats. Non fait : option `-U`, `-R`, sets d'amorces, mode B, sortie CSV, CLI `obiprimers`. + +## 14. Commande `obiprimers`, mode A (2026-11-03) + +Code : `pkg/obiprimer/search.go` (`SearchParams`, `DefaultSearchParams`, `Search`, `SearchInfo` ; +`FindPairs` renvoie maintenant aussi un `PairInfo`), `pkg/obitools/obiprimers/{options,obiprimers,output}.go`, +`cmd/obitools/obiprimers/main.go`, test `obiprimers_test.go`. Le `Makefile` découvre la commande tout seul +(`wildcard cmd/obitools/*/*.go`). Pattern : `InputOptionSet` + `OutputModeOptionSet` + `-t` obligatoire, +sortie par `obiitercsv` + `obicsv.CLICSVWriter` (donc `-o`, `--compress`), comme `obitaxonomy`. + +Options (noms alignés sur `obipcr` et `obigrep` quand ils existent) : + +| Option | ecoPrimers | Défaut | +|---|---|---| +| `-t/--taxonomy` (obligatoire) | `-d` (base) | | +| `-r/--restrict-to-taxon` | `-r` | tous les exemples | +| `-i/--ignore-taxon` | `-i` | | +| `-E/--exception-taxon` | `-E` | | +| `--rank` | `-t` | species | +| `-c/--circular` | `-c` | non | +| `--primer-length` | `-O` | 18 | +| `-e/--allowed-mismatches` | `-e` | 3 | +| `-l/--min-length`, `-L/--max-length` | `-l`, `-L` | 0, 1000 | +| `-3/--three-prime` | `-3` | 0 | +| `-q/--strict-quorum` | `-q` | 0,7 | +| `-s/--sensitivity` | `-s` | 0,9 | +| `-x/--false-positive` | `-x` | 0,1 | +| `-T/--specificity` | `-T` | 0,6 | +| `--salt-correction`, `--salt`, `--primer-concentration` | `-m`, `-a` | santalucia, 0,05 M, 8e-7 M | +| `--partitions` | (nouveau) | 256 | + +`-O` n'a pas d'alias court : dans obitools4 il désigne l'en-tête OBI en sortie. Non reprises : `-U` +(pas de multi-match), `-R` (séquence de référence), `-p/-b/-M/-k` (sets et liens), `-f`, `-v`, `-A`. + +Sortie CSV (en-têtes nommés, ordre des colonnes d'ecoPrimers) : `index, primer1, primer2, tm1, tm1_min, +tm2, tm2_min, gc1, gc2, good1, good2, examples, counterexamples, yule, example_taxa, +counterexample_taxa, bc, identified_taxa, bs, min_length, max_length, mean_length`. Tm arrondies à +0,1 °C, indices à 4 décimales, `NA` quand une Tm est indisponible. Le résumé (nombres de mots, de paires, +d'exemples...) va dans le journal (stderr), pas dans le CSV : les lignes `#` d'ecoPrimers sont donc +remplacées par des messages `INFO`. + +**Contrôle** : `obiprimers -t ncbitaxo_20250211.tgz -r 7953 -e 3 --three-prime 2 -l 30 -L 150 -c +mito.vert.fasta.gz` : 3 448 lignes, comparées à `Cyprinidae.ecoprimers` champ par champ (script jetable) : +**0 différence**. Teleostei : 407 lignes (les 6 écarts d'une séquence de la section 12 s'appliquent). + +Mesures du binaire complet (M3 Max, 16 workers, taxonomie NCBI complète comprise) : + +| Jeu | Temps | Pic mémoire | ecoPrimers (temps, pic) | +|---|---|---|---| +| Cyprinidae | 7,4 s | 2,82 Go | 134 s, 1,86 Go | +| Teleostei | 8,4 s | 3,42 Go | 350 s, 2,11 Go | + +**Le pic de mémoire est PLUS ÉLEVÉ qu'ecoPrimers**, malgré un temps 18 à 40 fois meilleur. Éléments mesurés : +le chargeur de taxonomie NCBI existant coûte à lui seul 1,07 Go (`obitaxonomy`, mesuré) ; avec +`GOGC=25` le pic de Cyprinidae passe de 2,82 à 1,98 Go (donc ~0,85 Go de ramasse-miettes en retard, au +prix de 2x de temps CPU) ; le reste (~0,9 Go au-dessus de la taxonomie pour 5,4 Mb d'exemples) n'est pas +analysé. Pistes non explorées : ne pas parser les 131 Mb du FASTA des séquences ignorées, libérer la +taxonomie après `AssignTaxa`, régler le GC de la commande. + +Avertissements `Taxid ... has to be updated` : viennent du mécanisme existant d'obitools (taxids +fusionnés dans le FASTA) ; `--update-taxid` les corrige à la lecture. + +Reste à faire : mode B, `-U`, `-R`, tests shell `obitests/obitools/obiprimers` (données dans le dépôt), +documentation, réduction de la mémoire, comparaison systématique `--rank` autre que species. diff --git a/cmd/obitools/obiprimers/main.go b/cmd/obitools/obiprimers/main.go new file mode 100644 index 0000000..29f095e --- /dev/null +++ b/cmd/obitools/obiprimers/main.go @@ -0,0 +1,41 @@ +package main + +import ( + "os" + + log "github.com/sirupsen/logrus" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obioptions" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obiconvert" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obicsv" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obiprimers" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiutils" +) + +func main() { + optionParser := obioptions.GenerateOptionParser( + "obiprimers", + "designs the primer pairs amplifying a set of example sequences and evaluates the DNA barcode they delimit", + obiprimers.OptionSet) + + _, args := optionParser(os.Args) + + if !obidefault.HasSelectedTaxonomy() { + log.Fatal("you must indicate a taxonomy using the -t or --taxonomy option") + } + + sequences, err := obiconvert.CLIReadBioSequences(args...) + if err != nil { + log.Fatalf("Cannot open the sequence files: %v", err) + } + + pairs, err := obiprimers.CLIPrimerPairs(sequences) + if err != nil { + log.Fatalf("%v", err) + } + + obicsv.CLICSVWriter(obiprimers.PairsToCSVIterator(pairs), true) + + obiutils.WaitForLastPipe() +} diff --git a/pkg/obiprimer/pairs_eval.go b/pkg/obiprimer/pairs_eval.go index ece6e0a..cd37864 100644 --- a/pkg/obiprimer/pairs_eval.go +++ b/pkg/obiprimer/pairs_eval.go @@ -36,6 +36,14 @@ type PairStats struct { // Score is the ranking criterion of ecoPrimers. func (p PairStats) Score() float32 { return p.Yule * p.Bs } +// PairInfo summarizes the pair search, for reporting. +type PairInfo struct { + Primers int // conserved words kept as candidate primers + 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 @@ -78,21 +86,29 @@ func ampHash(raw []byte, from, to int, reverse bool) uint64 { // is sorted by decreasing Yule*Bs. // // model computes the melting temperatures; nil leaves them NaN. -func FindPairs(c *Corpus, taxa *Taxa, wm []WordMatches, p PairParams, model *obithermo.Model) ([]PairStats, error) { +func FindPairs(c *Corpus, taxa *Taxa, wm []WordMatches, p PairParams, model *obithermo.Model) ([]PairStats, PairInfo, error) { + var info PairInfo if err := p.validate(); err != nil { - return nil, err + return nil, info, err } if taxa == nil || len(taxa.Of) != len(c.Seqs) { - return nil, fmt.Errorf("obiprimer: taxa do not match the corpus") + return nil, info, fmt.Errorf("obiprimer: taxa do not match the corpus") } ps := selectPrimers(c, wm, p) + info.Primers = len(ps.primers) + for _, pr := range ps.primers { + if pr.Good { + info.GoodPrimers++ + } + } if len(ps.primers) == 0 { - return nil, nil + return nil, info, nil } // Phase 1: which pairs amplify enough examples and few enough counterexamples. counts := c.Counts() acc := ps.countPairs(c, p) + info.CandidatePairs = len(acc) sens, fpq := float32(p.Sensitivity), float32(p.FalsePositive) var kept []pairKey for k, a := range acc { @@ -105,6 +121,7 @@ func FindPairs(c *Corpus, taxa *Taxa, wm []WordMatches, p PairParams, model *obi kept = append(kept, k) } } + info.RetainedPairs = len(kept) slices.SortFunc(kept, func(a, b pairKey) int { return cmp.Or(cmp.Compare(a.p1, b.p1), cmp.Compare(a.p2, b.p2), boolCmp(a.a1, b.a1), boolCmp(a.a2, b.a2)) @@ -131,7 +148,7 @@ func FindPairs(c *Corpus, taxa *Taxa, wm []WordMatches, p PairParams, model *obi slices.SortFunc(out, func(a, b PairStats) int { return cmp.Or(-cmp.Compare(a.Score(), b.Score()), cmp.Compare(a.Primer1, b.Primer1), cmp.Compare(a.Primer2, b.Primer2)) }) - return out, nil + return out, info, nil } func boolCmp(a, b bool) int { diff --git a/pkg/obiprimer/pairs_test.go b/pkg/obiprimer/pairs_test.go index e167772..9160cd7 100644 --- a/pkg/obiprimer/pairs_test.go +++ b/pkg/obiprimer/pairs_test.go @@ -141,7 +141,7 @@ func TestPairsAgainstEcoPrimers(t *testing.T) { t.Fatal(err) } lap("approximate search") - got, err := FindPairs(corpus, taxa, wm, PairParams{ + got, _, err := FindPairs(corpus, taxa, wm, PairParams{ PrimerLength: 18, MinLength: 30, MaxLength: 150, ThreePrime: 2, Sensitivity: 0.9, FalsePositive: 0.1, MinSpecificity: 0.6, Circular: true, }, obithermo.DefaultModel()) diff --git a/pkg/obiprimer/search.go b/pkg/obiprimer/search.go new file mode 100644 index 0000000..3a48567 --- /dev/null +++ b/pkg/obiprimer/search.go @@ -0,0 +1,91 @@ +package obiprimer + +import ( + "fmt" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obithermo" +) + +// SearchParams gathers every parameter of the primer pair search, named after +// the ecoPrimers options they come from. +type SearchParams struct { + PrimerLength int // -O: size of the primers (default 18) + MaxErrors int // -e: substitutions allowed between a primer and its sites (default 3) + MinLength int // -l: the amplicon (primers excluded) must be longer than this + MaxLength int // -L: and shorter than this + ThreePrime int // -3: number of 3' bases that must match exactly + StrictQuorum float64 // -q: fraction of the examples containing the exact word (0.7) + Sensitivity float64 // -s: fraction of the examples a primer must match (0.9) + FalsePositive float64 // -x: fraction of the counterexamples it may match (0.1) + MinSpecificity float64 // -T: pairs with a lower Bs are dropped (0.6) + Circular bool // -c: sequences are circular + + Partitions int // partitions of the k-mer counter (default 256) + TmpDir string + + Model *obithermo.Model // melting temperatures; nil leaves them NaN +} + +// DefaultSearchParams returns the ecoPrimers defaults. +func DefaultSearchParams() SearchParams { + return SearchParams{ + PrimerLength: 18, MaxErrors: 3, MinLength: 0, MaxLength: 1000, + StrictQuorum: 0.7, Sensitivity: 0.9, FalsePositive: 0.1, MinSpecificity: 0.6, + Partitions: 256, Model: obithermo.DefaultModel(), + } +} + +// SearchInfo reports the size of every stage of a search. +type SearchInfo struct { + Counts Counts + StrictWords int // words found without error in at least StrictQuorum of the examples + MinStrict int // that quorum, in sequences + Matches int // approximate occurrences of those words + Pair PairInfo +} + +// Search runs the whole ecoPrimers pipeline on a corpus: strict words, +// approximate occurrences, pairs and their indices. The pairs are sorted by +// decreasing Yule*Bs. +func Search(c *Corpus, taxa *Taxa, p SearchParams) ([]PairStats, *SearchInfo, error) { + if p.MaxErrors < 0 { + return nil, nil, fmt.Errorf("obiprimer: the allowed mismatch count cannot be negative") + } + + sw, err := FindStrictWords(c, WordSearch{ + PrimerLength: p.PrimerLength, + Quorum: p.StrictQuorum, + Circular: p.Circular, + Partitions: p.Partitions, + TmpDir: p.TmpDir, + }) + if err != nil { + return nil, nil, err + } + info := &SearchInfo{Counts: sw.Counts, StrictWords: len(sw.Words), MinStrict: sw.MinSequences} + + words := make([]uint64, len(sw.Words)) + for i, w := range sw.Words { + words[i] = w.Kmer + } + wm, err := MatchWords(c, words, MatchParams{K: p.PrimerLength, MaxErrors: p.MaxErrors, Circular: p.Circular}) + if err != nil { + return nil, info, err + } + for _, w := range wm { + info.Matches += len(w.Matches) + } + + pairs, pinfo, err := FindPairs(c, taxa, wm, PairParams{ + PrimerLength: p.PrimerLength, + MinLength: p.MinLength, + MaxLength: p.MaxLength, + ThreePrime: p.ThreePrime, + Sensitivity: p.Sensitivity, + FalsePositive: p.FalsePositive, + MinSpecificity: p.MinSpecificity, + Circular: p.Circular, + }, p.Model) + info.Pair = pinfo + return pairs, info, err +} diff --git a/pkg/obitools/obiprimers/obiprimers.go b/pkg/obitools/obiprimers/obiprimers.go new file mode 100644 index 0000000..0362434 --- /dev/null +++ b/pkg/obitools/obiprimers/obiprimers.go @@ -0,0 +1,62 @@ +package obiprimers + +import ( + "fmt" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiiter" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiprimer" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitax" + log "github.com/sirupsen/logrus" +) + +// CLIPrimerPairs reads the sequences, selects the examples and the +// counterexamples from the taxonomy options, and searches the primer pairs. +// The default taxonomy must be loaded. +func CLIPrimerPairs(sequences obiiter.IBioSequence) ([]obiprimer.PairStats, error) { + taxonomy := obitax.DefaultTaxonomy() + if taxonomy == nil { + return nil, fmt.Errorf("a taxonomy is required: use the -t option") + } + + sel := obiprimer.Selection{} + var err error + if sel.Restrict, err = obiprimer.BelongingTo(taxonomy, CLIRestrictTaxa()); err != nil { + return nil, err + } + if sel.Ignore, err = obiprimer.BelongingTo(taxonomy, CLIIgnoreTaxa()); err != nil { + return nil, err + } + if sel.Except, err = obiprimer.BelongingTo(taxonomy, CLIExceptionTaxa()); err != nil { + return nil, err + } + + params, err := CLISearchParams() + if err != nil { + return nil, err + } + + log.Infof("Reading the sequences") + corpus, err := obiprimer.LoadCorpus(sequences, sel) + if err != nil { + return nil, err + } + taxa, err := corpus.AssignTaxa(taxonomy, CLIRank()) + if err != nil { + return nil, err + } + + counts := corpus.Counts() + log.Infof("%d example sequences (%d %s), %d counterexamples (%d %s), %d ignored", + counts.Examples, taxa.ExampleTaxa, CLIRank(), + counts.CounterExamples, taxa.CounterExamples, CLIRank(), counts.Ignored) + + pairs, info, err := obiprimer.Search(corpus, taxa, params) + if info != nil { + log.Infof("%d strict words (found without error in at least %d examples)", info.StrictWords, info.MinStrict) + 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)) + } + return pairs, err +} diff --git a/pkg/obitools/obiprimers/obiprimers_test.go b/pkg/obitools/obiprimers/obiprimers_test.go new file mode 100644 index 0000000..9c33a27 --- /dev/null +++ b/pkg/obitools/obiprimers/obiprimers_test.go @@ -0,0 +1,115 @@ +package obiprimers + +import ( + "math" + "testing" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiprimer" +) + +func samplePairs() []obiprimer.PairStats { + return []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, + 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()}, + } +} + +func TestPairsToCSVIterator(t *testing.T) { + it := PairsToCSVIterator(samplePairs()) + + header := it.Header() + if len(header) != len(csvColumns) || header[0] != "index" || header[1] != "primer1" { + t.Fatalf("unexpected header %v", header) + } + + var records []map[string]interface{} + for it.Next() { + batch := it.Get() + for _, r := range batch.Slice() { + records = append(records, r) + } + } + if len(records) != 2 { + t.Fatalf("%d records, want 2", len(records)) + } + + r := records[0] + 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["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"]) + } + if _, present := r["tm2_min"]; present { + t.Error("an unavailable temperature must be left out, to be written as NA") + } + if _, present := records[1]["tm1"]; present { + t.Error("NaN temperatures must be left out") + } + if records[1]["index"] != 2 { + t.Errorf("second pair index %v", records[1]["index"]) + } +} + +func TestEmptyResultStillHasAHeader(t *testing.T) { + it := PairsToCSVIterator(nil) + if len(it.Header()) != len(csvColumns) { + t.Fatalf("header %v", it.Header()) + } + n := 0 + for it.Next() { + batch := it.Get() + n += len(batch.Slice()) + } + if n != 0 { + t.Errorf("%d records for no pairs", n) + } +} + +func TestEveryColumnIsProduced(t *testing.T) { + full := pairRecord(1, samplePairs()[0]) + full["tm2_min"] = 1.0 // the only column left out for this pair + for _, c := range csvColumns { + if _, ok := full[c]; !ok { + t.Errorf("column %q is never filled", c) + } + } + if len(full) != len(csvColumns) { + t.Errorf("record has %d fields for %d columns", len(full), len(csvColumns)) + } +} + +func TestCLISearchParams(t *testing.T) { + p, err := CLISearchParams() + if err != nil { + t.Fatal(err) + } + d := obiprimer.DefaultSearchParams() + if p.PrimerLength != d.PrimerLength || p.MaxErrors != d.MaxErrors || p.MaxLength != d.MaxLength || + p.StrictQuorum != d.StrictQuorum || p.Sensitivity != d.Sensitivity || + p.FalsePositive != d.FalsePositive || p.MinSpecificity != d.MinSpecificity || p.Model == nil { + t.Errorf("option defaults %+v differ from the ecoPrimers defaults %+v", p, d) + } + + _saltMethod = "owczarzy" + if p, err = CLISearchParams(); err != nil || p.Model == nil { + t.Errorf("owczarzy must be accepted: %v", err) + } + _saltMethod = "nope" + if _, err = CLISearchParams(); err == nil { + t.Error("an unknown salt correction must be refused") + } + _saltMethod = "santalucia" + _salt = -1 + if _, err = CLISearchParams(); err == nil { + t.Error("a negative salt concentration must be refused") + } +} diff --git a/pkg/obitools/obiprimers/options.go b/pkg/obitools/obiprimers/options.go new file mode 100644 index 0000000..84cb611 --- /dev/null +++ b/pkg/obitools/obiprimers/options.go @@ -0,0 +1,182 @@ +// obiprimers function utility package. +// +// The obitools/obiprimers package contains every function specific to the +// obiprimers utility, the OBITools4 version of ecoPrimers: it designs the +// primer pairs that best amplify a set of example sequences, optionally +// avoiding a set of counterexample sequences, and measures the quality of the +// DNA marker each pair delimits. +package obiprimers + +import ( + "fmt" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obioptions" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiprimer" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obithermo" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obiconvert" + "github.com/DavidGamba/go-getoptions" +) + +var _restrict = []string{} +var _ignore = []string{} +var _except = []string{} +var _rank = "species" + +var _primerLength = 18 +var _allowedMismatch = 3 +var _minLength = 0 +var _maxLength = 1000 +var _threePrime = 0 +var _circular = false + +var _strictQuorum = 0.7 +var _sensitivity = 0.9 +var _falsePositive = 0.1 +var _specificity = 0.6 + +var _saltMethod = "santalucia" +var _salt = obithermo.DefaultSalt +var _primerConc = obithermo.DefaultPrimerConc + +var _partitions = 256 + +// OptionSet registers the options of obiprimers. +func OptionSet(options *getoptions.GetOpt) { + obiconvert.InputOptionSet(options) + obiconvert.OutputModeOptionSet(options, true) + obioptions.LoadTaxonomyOptionSet(options, true, false) + + 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.StringSliceVar(&_ignore, "ignore-taxon", 1, 1, + options.Alias("i"), + options.ArgName("TAXID"), + options.Description("Declares the sequences belonging to this taxon as counterexamples "+ + "(sequences the primers must not amplify). Can be repeated.")) + + options.StringSliceVar(&_except, "exception-taxon", 1, 1, + options.Alias("E"), + options.ArgName("TAXID"), + options.Description("Declares a sub-clade of the examples as counterexamples. Can be repeated.")) + + options.StringVar(&_rank, "rank", _rank, + options.ArgName("RANK"), + options.Description("Taxonomic rank at which coverage (Bc) and specificity (Bs) are computed.")) + + options.BoolVar(&_circular, "circular", _circular, + options.Alias("c"), + options.Description("Considers the sequences as circular (whole organelle genomes).")) + + options.IntVar(&_primerLength, "primer-length", _primerLength, + options.ArgName("LENGTH"), + options.Description("Size of the primers, from 2 to 31.")) + + options.IntVar(&_allowedMismatch, "allowed-mismatches", _allowedMismatch, + options.Alias("e"), + options.ArgName("ERRORS"), + options.Description("Maximum number of substitutions between a primer and a site it matches. "+ + "The search is fast up to about 5; beyond, the number of matches explodes.")) + + options.IntVar(&_minLength, "min-length", _minLength, + options.Alias("l"), + options.ArgName("LENGTH"), + options.Description("The amplified region, primers excluded, must be longer than this.")) + + options.IntVar(&_maxLength, "max-length", _maxLength, + options.Alias("L"), + options.ArgName("LENGTH"), + options.Description("The amplified region, primers excluded, must be shorter than this.")) + + options.IntVar(&_threePrime, "three-prime", _threePrime, + options.Alias("3"), + options.ArgName("BASES"), + options.Description("Number of bases at the 3' end of each primer that must match the "+ + "amplified sequence exactly (an ambiguous base is not a match).")) + + options.Float64Var(&_strictQuorum, "strict-quorum", _strictQuorum, + options.Alias("q"), + options.ArgName("FRACTION"), + options.Description("Fraction of the examples in which a primer word must be found without any error.")) + + options.Float64Var(&_sensitivity, "sensitivity", _sensitivity, + options.Alias("s"), + options.ArgName("FRACTION"), + options.Description("Fraction of the examples a primer must match (with allowed mismatches), "+ + "and a pair must amplify.")) + + options.Float64Var(&_falsePositive, "false-positive", _falsePositive, + options.Alias("x"), + options.ArgName("FRACTION"), + options.Description("Fraction of the counterexamples a good primer may match, and a pair may amplify.")) + + options.Float64Var(&_specificity, "specificity", _specificity, + options.Alias("T"), + options.ArgName("FRACTION"), + options.Description("Pairs whose specificity Bs is lower than this are dropped.")) + + options.StringVar(&_saltMethod, "salt-correction", _saltMethod, + options.ArgName("METHOD"), + options.Description("Salt correction of the melting temperatures: 'santalucia' or 'owczarzy'.")) + + options.Float64Var(&_salt, "salt", _salt, + options.ArgName("MOLAR"), + options.Description("Monovalent salt concentration (M) used for the melting temperatures.")) + + options.Float64Var(&_primerConc, "primer-concentration", _primerConc, + options.ArgName("MOLAR"), + options.Description("Primer concentration (M) used for the melting temperatures.")) + + options.IntVar(&_partitions, "partitions", _partitions, + options.ArgName("COUNT"), + options.Description("Number of temporary partition files used to count the conserved words. "+ + "Increase it for very large inputs (memory per partition grows with data size / COUNT).")) +} + +// CLIRestrictTaxa returns the taxids restricting the examples. +func CLIRestrictTaxa() []string { return _restrict } + +// CLIIgnoreTaxa returns the taxids of the counterexample clades. +func CLIIgnoreTaxa() []string { return _ignore } + +// CLIExceptionTaxa returns the taxids of the sub-clades turned into counterexamples. +func CLIExceptionTaxa() []string { return _except } + +// CLIRank returns the taxonomic rank of the coverage and specificity indices. +func CLIRank() string { return _rank } + +// CLISearchParams builds the search parameters from the command line. +func CLISearchParams() (obiprimer.SearchParams, error) { + var method obithermo.SaltMethod + switch _saltMethod { + case "santalucia": + method = obithermo.SantaLucia + case "owczarzy": + method = obithermo.Owczarzy + default: + return obiprimer.SearchParams{}, fmt.Errorf("unknown salt correction %q: use 'santalucia' or 'owczarzy'", _saltMethod) + } + model, err := obithermo.NewModel(_primerConc, obithermo.DefaultTargetConc, _salt, method) + if err != nil { + return obiprimer.SearchParams{}, err + } + + return obiprimer.SearchParams{ + PrimerLength: _primerLength, + MaxErrors: _allowedMismatch, + MinLength: _minLength, + MaxLength: _maxLength, + ThreePrime: _threePrime, + StrictQuorum: _strictQuorum, + Sensitivity: _sensitivity, + FalsePositive: _falsePositive, + MinSpecificity: _specificity, + Circular: _circular, + Partitions: _partitions, + Model: model, + }, nil +} diff --git a/pkg/obitools/obiprimers/output.go b/pkg/obitools/obiprimers/output.go new file mode 100644 index 0000000..49c3dc6 --- /dev/null +++ b/pkg/obitools/obiprimers/output.go @@ -0,0 +1,78 @@ +package obiprimers + +import ( + "math" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiitercsv" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiprimer" +) + +// Columns of the output, in the order of the ecoPrimers result table. +var csvColumns = []string{ + "index", // 1 serial number + "primer1", "primer2", // 2-3 + "tm1", "tm1_min", "tm2", "tm2_min", // 4-7 melting temperatures: perfect match, lowest over the examples + "gc1", "gc2", // 8-9 + "good1", "good2", // 10 primer conserved in the examples and rare in the counterexamples + "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 + "min_length", "max_length", "mean_length", // 19-21 +} + +// round returns x rounded to the given number of decimals; the CSV writer +// prints it with its shortest representation. +func round(x float64, decimals int) float64 { + p := math.Pow10(decimals) + return math.Round(x*p) / p +} + +// pairRecord converts a pair to a CSV record. Unavailable values (NaN) are +// left out, which the writer prints as the NA value. +func pairRecord(index int, p obiprimer.PairStats) obiitercsv.CSVRecord { + r := obiitercsv.CSVRecord{ + "index": index, "primer1": p.Primer1, "primer2": p.Primer2, + "gc1": p.GC1, "gc2": p.GC2, "good1": p.Good1, "good2": p.Good2, + "examples": p.Examples, "counterexamples": p.CounterExamples, + "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), + "min_length": p.MinLen, "max_length": p.MaxLen, "mean_length": round(p.MeanLen, 2), + } + for key, v := range map[string]float64{"tm1": p.Tm1, "tm1_min": p.MinTm1, "tm2": p.Tm2, "tm2_min": p.MinTm2} { + if !math.IsNaN(v) { + r[key] = round(v, 1) + } + } + return r +} + +// PairsToCSVIterator streams the pairs, in order, as CSV records. +func PairsToCSVIterator(pairs []obiprimer.PairStats) *obiitercsv.ICSVRecord { + iterator := obiitercsv.NewICSVRecord() + iterator.Add(1) + for _, column := range csvColumns { + iterator.AppendField(column) + } + + go func() { + defer iterator.Done() + batchSize := obidefault.BatchSize() + order := 0 + for start := 0; start < len(pairs) || order == 0; start += batchSize { + end := min(start+batchSize, len(pairs)) + data := make([]obiitercsv.CSVRecord, 0, end-start) + for i := start; i < end; i++ { + data = append(data, pairRecord(i+1, pairs[i])) + } + iterator.Push(obiitercsv.MakeCSVRecordBatch("obiprimers", order, data)) + order++ + } + iterator.Close() + }() + return iterator +}