From 3087c70daa1e2c33172e923a4e5aa00a52ae7d43 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Mon, 2 Nov 2026 22:12:43 +1100 Subject: [PATCH] Feature: Implement constrained k-mer generation and clique search This commit introduces core functionality for oligonucleotide family tagging, constrained k-mer generation, and clique finding. It also includes updates to the Go language version and refactors the logic for k-mer dissimilarity metrics. --- Release-notes.md | 17 ++ cmd/obitools/obioligotag/main.go | 33 +++ go.mod | 2 +- go.work | 2 +- obitests/obitools/obioligotag/test.sh | 205 ++++++++++++++++++ pkg/obigraph/clique.go | 176 +++++++++++++++ pkg/obigraph/clique_test.go | 143 ++++++++++++ pkg/obikmer/enumkmer.go | 60 ++++- pkg/obikmer/kmer_constraints.go | 141 ++++++++++++ pkg/obikmer/kmer_constraints_test.go | 123 +++++++++++ pkg/obikmer/kmer_dissimilarity.go | 84 +++++++ pkg/obikmer/kmer_dissimilarity_test.go | 113 ++++++++++ pkg/obiseq/iupac_nog.go | 13 ++ pkg/obitools/obioligotag/input.go | 148 +++++++++++++ pkg/obitools/obioligotag/obioligotag.go | 217 +++++++++++++++++++ pkg/obitools/obioligotag/obioligotag_test.go | 177 +++++++++++++++ pkg/obitools/obioligotag/options.go | 150 +++++++++++++ 17 files changed, 1791 insertions(+), 13 deletions(-) create mode 100644 cmd/obitools/obioligotag/main.go create mode 100755 obitests/obitools/obioligotag/test.sh create mode 100644 pkg/obigraph/clique.go create mode 100644 pkg/obigraph/clique_test.go create mode 100644 pkg/obikmer/kmer_constraints.go create mode 100644 pkg/obikmer/kmer_constraints_test.go create mode 100644 pkg/obikmer/kmer_dissimilarity.go create mode 100644 pkg/obikmer/kmer_dissimilarity_test.go create mode 100644 pkg/obitools/obioligotag/input.go create mode 100644 pkg/obitools/obioligotag/obioligotag.go create mode 100644 pkg/obitools/obioligotag/obioligotag_test.go create mode 100644 pkg/obitools/obioligotag/options.go diff --git a/Release-notes.md b/Release-notes.md index 143576a..1666e24 100644 --- a/Release-notes.md +++ b/Release-notes.md @@ -2,6 +2,23 @@ ## New changes +### New features + +- New `obioligotag` command, a re-implementation of the `oligotag` command of + OBITools 2. It builds a family of oligonucleotides (tags) of at least + `--family-size` members that are pairwise at a distance of at least + `--distance`. Candidates are either generated (`--oligo-size`, up to 12) or + read from a file (`--oligo-list`, sequence file or one oligonucleotide per + line). They can be restricted by alphabet (`--alphabet`, for instance `ACT` + to forbid G), GC content, homopolymers and IUPAC patterns. The distance is + either the Hamming distance or the LCS-mode alignment distance + (`--distance-type`). + +- New functions in `obikmer` to compare and enumerate k-mers encoded as + `uint64`: `HammingDistance`, `KmerLCS` and `KmerLCSBatch` (SIMD when built with + `GOEXPERIMENT=simd`), `IterConstrainedKmers`, and in `obigraph` a clique search + (`FindClique`) that does not need the adjacency matrix. + ### Bug fixes - In `obipairing` correct the misspelling of the `obiparing_*` tags where the `i` diff --git a/cmd/obitools/obioligotag/main.go b/cmd/obitools/obioligotag/main.go new file mode 100644 index 0000000..04cce0f --- /dev/null +++ b/cmd/obitools/obioligotag/main.go @@ -0,0 +1,33 @@ +package main + +import ( + "os" + + log "github.com/sirupsen/logrus" + + "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/obioligotag" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiutils" +) + +func main() { + optionParser := obioptions.GenerateOptionParser( + "obioligotag", + "designs a family of oligonucleotides (tags) at a minimum distance from each other", + obioligotag.OptionSet) + + _, args := optionParser(os.Args) + if len(args) > 0 { + log.Fatalf("unexpected arguments %v: use --oligo-list to provide candidate oligonucleotides", args) + } + + oligos, err := obioligotag.CLIOligoFamily() + if err != nil { + log.Fatalf("%v", err) + } + + obiconvert.CLIWriteBioSequences(obioligotag.CLIOligosToIterator(oligos), true) + + obiutils.WaitForLastPipe() +} diff --git a/go.mod b/go.mod index 86a7baf..36dd212 100644 --- a/go.mod +++ b/go.mod @@ -1,6 +1,6 @@ module git.metabarcoding.org/obitools/obitools4/obitools4 -go 1.26.1 +go 1.27 require ( github.com/DavidGamba/go-getoptions v0.33.0 diff --git a/go.work b/go.work index e4243b4..47208ec 100644 --- a/go.work +++ b/go.work @@ -1,3 +1,3 @@ -go 1.26.1 +go 1.27 use . diff --git a/obitests/obitools/obioligotag/test.sh b/obitests/obitools/obioligotag/test.sh new file mode 100755 index 0000000..86fd52b --- /dev/null +++ b/obitests/obitools/obioligotag/test.sh @@ -0,0 +1,205 @@ +#!/bin/bash + +# +# Here give the name of the test serie +# + +TEST_NAME=obioligotag +CMD=obioligotag + +###### +# +# Some variable and function definitions: please don't change them +# +###### +TEST_DIR="$(dirname "$(readlink -f "${BASH_SOURCE[0]}")")" +OBITOOLS_DIR="${TEST_DIR/obitest*/}build" +export PATH="${OBITOOLS_DIR}:${PATH}" + +MCMD="$(echo "${CMD:0:4}" | tr '[:lower:]' '[:upper:]')$(echo "${CMD:4}" | tr '[:upper:]' '[:lower:]')" + +TMPDIR="$(mktemp -d)" +ntest=0 +success=0 +failed=0 + +cleanup() { + echo "========================================" 1>&2 + echo "## Results of the $TEST_NAME tests:" 1>&2 + + echo 1>&2 + echo "- $ntest tests run" 1>&2 + echo "- $success successfully completed" 1>&2 + echo "- $failed failed tests" 1>&2 + echo 1>&2 + echo "Cleaning up the temporary directory..." 1>&2 + echo 1>&2 + echo "========================================" 1>&2 + + rm -rf "$TMPDIR" # Suppress the temporary directory + + if [ $failed -gt 0 ]; then + log "$TEST_NAME tests failed" + log + log + exit 1 + fi + + log + log + + exit 0 +} + +log() { + echo -e "[$TEST_NAME @ $(date)] $*" 1>&2 +} + +log "Testing $TEST_NAME..." +log "Test directory is $TEST_DIR" +log "obitools directory is $OBITOOLS_DIR" +log "Temporary directory is $TMPDIR" +log "files: $(find $TEST_DIR | awk -F'/' '{print $NF}' | tail -n +2)" + +###################################################################### +#### +#### Below are the tests +#### +#### Before each test : +#### - increment the variable ntest +#### +#### Run the command as the condition of an if / then /else +#### - The command must return 0 on success +#### - The command must return an exit code different from 0 on failure +#### - The datafiles are stored in the same directory than the test script +#### - The test script directory is stored in the TEST_DIR variable +#### - If result files have to be produced they must be stored +#### in the temporary directory (TMPDIR variable) +#### +#### then clause is executed on success of the command +#### - Write a success message using the log function +#### - increment the variable success +#### +#### else clause is executed on failure of the command +#### - Write a failure message using the log function +#### - increment the variable failed +#### +###################################################################### + + + +# +# Every result is checked on its properties (size, alphabet, constraints), +# not against a reference file, because several equivalent families exist. +# + +printf 'acgtac\nacgtat\nttttgg\nggcaac\ncaccat\nacgtac\n' > "${TMPDIR}/oligos.txt" +printf '>a\nACGTAC\n>b\nACGTAT\n>c\nTTTTGG\n>d\nGGCAAC\n>e\nCACCAT\n' > "${TMPDIR}/oligos.fasta" + +((ntest++)) +if $CMD -h > "${TMPDIR}/help.txt" 2>&1 +then + log "$MCMD: printing help OK" + ((success++)) +else + log "$MCMD: printing help failed" + ((failed++)) +fi + +((ntest++)) +if $CMD -s 5 -f 24 -d 3 -g 3 -p 2 > "${TMPDIR}/tags.fasta" 2> /dev/null \ + && [ "$(grep -vc '>' "${TMPDIR}/tags.fasta")" -ge 24 ] \ + && ! grep -v '>' "${TMPDIR}/tags.fasta" | grep -qE '(a{3}|c{3}|g{3}|t{3})' +then + log "$MCMD: generating a family with homopolymer constraint OK" + ((success++)) +else + log "$MCMD: generating a family with homopolymer constraint failed" + ((failed++)) +fi + +((ntest++)) +if $CMD -s 7 -f 60 -d 3 --alphabet ACT > "${TMPDIR}/noG.fasta" 2> /dev/null \ + && [ "$(grep -vc '>' "${TMPDIR}/noG.fasta")" -ge 60 ] \ + && ! grep -v '>' "${TMPDIR}/noG.fasta" | grep -q 'g' +then + log "$MCMD: restricted alphabet OK" + ((success++)) +else + log "$MCMD: restricted alphabet failed" + ((failed++)) +fi + +((ntest++)) +if $CMD -s 6 -f 20 -d 3 --distance-type lcs > "${TMPDIR}/lcs.fasta" 2> /dev/null \ + && [ "$(grep -vc '>' "${TMPDIR}/lcs.fasta")" -ge 20 ] +then + log "$MCMD: LCS distance OK" + ((success++)) +else + log "$MCMD: LCS distance failed" + ((failed++)) +fi + +((ntest++)) +if $CMD -s 7 -f 40 -d 3 -r cnnnnnn > "${TMPDIR}/rejected.fasta" 2> /dev/null \ + && [ "$(grep -vc '>' "${TMPDIR}/rejected.fasta")" -ge 40 ] \ + && ! grep -v '>' "${TMPDIR}/rejected.fasta" | grep -q '^c' +then + log "$MCMD: rejected pattern OK" + ((success++)) +else + log "$MCMD: rejected pattern failed" + ((failed++)) +fi + +((ntest++)) +if $CMD -L "${TMPDIR}/oligos.txt" -f 3 -d 3 > "${TMPDIR}/from_text.fasta" 2> /dev/null \ + && [ "$(grep -vc '>' "${TMPDIR}/from_text.fasta")" -ge 3 ] +then + log "$MCMD: selecting oligonucleotides from a text file OK" + ((success++)) +else + log "$MCMD: selecting oligonucleotides from a text file failed" + ((failed++)) +fi + +((ntest++)) +if $CMD -L "${TMPDIR}/oligos.fasta" -f 3 -d 3 > "${TMPDIR}/from_fasta.fasta" 2> /dev/null \ + && [ "$(grep -vc '>' "${TMPDIR}/from_fasta.fasta")" -ge 3 ] +then + log "$MCMD: selecting oligonucleotides from a fasta file OK" + ((success++)) +else + log "$MCMD: selecting oligonucleotides from a fasta file failed" + ((failed++)) +fi + +((ntest++)) +if ! $CMD -s 5 -f 100 -d 3 > /dev/null 2>&1 +then + log "$MCMD: family size beyond the Hamming bound refused OK" + ((success++)) +else + log "$MCMD: family size beyond the Hamming bound not refused" + ((failed++)) +fi + +((ntest++)) +if ! $CMD -s 5 -L "${TMPDIR}/oligos.txt" -f 3 > /dev/null 2>&1 +then + log "$MCMD: --oligo-size with --oligo-list refused OK" + ((success++)) +else + log "$MCMD: --oligo-size with --oligo-list not refused" + ((failed++)) +fi + +######################################### +# +# At the end of the tests +# the cleanup function is called +# +######################################### + +cleanup diff --git a/pkg/obigraph/clique.go b/pkg/obigraph/clique.go new file mode 100644 index 0000000..8572dc0 --- /dev/null +++ b/pkg/obigraph/clique.go @@ -0,0 +1,176 @@ +package obigraph + +import ( + "context" + "math/rand" + "sync" +) + +// NeighborFilter defines an implicit graph over the vertices 0..n-1. +// It appends to out[:0] the vertices of candidates that are adjacent to v, +// keeping their order, and returns the resulting slice. It is called +// concurrently by several workers and must not modify candidates. +// +// Describing the adjacency by a filter, instead of a matrix, keeps the memory +// in O(n) and lets the caller evaluate adjacency in bulk (SIMD, for instance). +type NeighborFilter func(v int32, candidates, out []int32) []int32 + +// CliqueOptions tunes FindClique. +type CliqueOptions struct { + // Workers is the number of concurrent searches (default 1). + Workers int + // Seed seeds the random vertex orders of the restarts. + Seed int64 +} + +// CliqueResult is the outcome of FindClique. +type CliqueResult struct { + // Clique is the largest clique found (a clique of size >= target if Reached). + Clique []int32 + // Reached is true if len(Clique) >= target. + Reached bool + // Exhausted is true if the whole search space has been explored. If Reached + // is false, it proves that no clique of the target size exists. + Exhausted bool +} + +// FindClique looks for a clique of at least target vertices in the graph of n +// vertices described by filter. It stops as soon as one is found, when the +// search space is exhausted, or when ctx is done; in the last case, the largest +// clique met so far is returned. +// +// The search is a depth first search with backtracking, pruned when the +// current clique plus the remaining candidates cannot reach target. The first +// pass of the first worker follows the natural order of the vertices (greedy +// "lexicographic" construction); every other pass restarts from a random +// order, with a work budget that doubles at each restart, so that the search +// is complete in the long run. +func FindClique(ctx context.Context, n int, filter NeighborFilter, target int, opts CliqueOptions) CliqueResult { + if target <= 0 { + return CliqueResult{Reached: true, Exhausted: true} + } + if n <= 0 { + return CliqueResult{Exhausted: true} + } + + workers := max(opts.Workers, 1) + ctx, cancel := context.WithCancel(ctx) + defer cancel() + + var ( + mutex sync.Mutex + best []int32 + exhausted bool + wg sync.WaitGroup + ) + + publish := func(clique []int32) { + mutex.Lock() + if len(clique) > len(best) { + best = append(best[:0], clique...) + } + mutex.Unlock() + } + + for w := range workers { + wg.Go(func() { + s := &cliqueSearch{ + ctx: ctx, + filter: filter, + target: target, + publish: func(c []int32) { + publish(c) + if len(c) >= target { + cancel() + } + }, + } + rng := rand.New(rand.NewSource(opts.Seed + int64(w))) + candidates := make([]int32, n) + for i := range candidates { + candidates[i] = int32(i) + } + + for restart := 0; ctx.Err() == nil; restart++ { + if w > 0 || restart > 0 { + rng.Shuffle(n, func(i, j int) { candidates[i], candidates[j] = candidates[j], candidates[i] }) + } + s.budget = int64(n) * int64(target) << min(restart, 20) + s.work = 0 + s.aborted = false + s.search(0, candidates) + + if !s.aborted && ctx.Err() == nil { + mutex.Lock() + exhausted = true + mutex.Unlock() + cancel() + } + } + }) + } + wg.Wait() + + return CliqueResult{ + Clique: best, + Reached: len(best) >= target, + Exhausted: exhausted, + } +} + +type cliqueSearch struct { + ctx context.Context + filter NeighborFilter + target int + publish func([]int32) + + current []int32 // the clique under construction + buffers [][]int32 // candidate list of each depth + bestLen int + + budget int64 // maximum work (number of candidates tested) for this restart + work int64 + aborted bool // budget or context exhausted before the end of the exploration +} + +func (s *cliqueSearch) search(depth int, candidates []int32) { + for i := range candidates { + if depth+len(candidates)-i < s.target { + return + } + if s.work > s.budget || s.ctx.Err() != nil { + s.aborted = true + return + } + + v := candidates[i] + s.current = append(s.current[:depth], v) + if len(s.current) > s.bestLen { + s.bestLen = len(s.current) + s.publish(s.current) + if s.bestLen >= s.target { + s.aborted = true // the search is over: do not report an exhaustion + return + } + } + + rest := candidates[i+1:] + if depth+1+len(rest) < s.target { + continue + } + for len(s.buffers) <= depth+1 { + s.buffers = append(s.buffers, nil) + } + next := s.filter(v, rest, s.buffers[depth+1][:0]) + s.buffers[depth+1] = next + s.work += int64(len(rest)) + + if depth+1+len(next) < s.target { + continue + } + s.search(depth+1, next) + if s.aborted { + return + } + } +} diff --git a/pkg/obigraph/clique_test.go b/pkg/obigraph/clique_test.go new file mode 100644 index 0000000..b2c6cd2 --- /dev/null +++ b/pkg/obigraph/clique_test.go @@ -0,0 +1,143 @@ +package obigraph + +import ( + "context" + "math/bits" + "math/rand" + "testing" + "time" +) + +// adjacency matrix filter for small graphs +func matrixFilter(adj [][]bool) NeighborFilter { + return func(v int32, candidates, out []int32) []int32 { + for _, c := range candidates { + if adj[v][c] { + out = append(out, c) + } + } + return out + } +} + +func randomGraph(rng *rand.Rand, n int, p float64) [][]bool { + adj := make([][]bool, n) + for i := range adj { + adj[i] = make([]bool, n) + } + for i := 0; i < n; i++ { + for j := i + 1; j < n; j++ { + if rng.Float64() < p { + adj[i][j], adj[j][i] = true, true + } + } + } + return adj +} + +// maximum clique size by brute force over all the subsets (n <= 20) +func bruteMaxClique(adj [][]bool) int { + n := len(adj) + best := 0 + for mask := 1; mask < 1<>i&1 == 0 { + continue + } + for j := i + 1; j < n; j++ { + if mask>>j&1 == 1 && !adj[i][j] { + ok = false + break + } + } + } + if ok { + best = size + } + } + return best +} + +func isClique(adj [][]bool, c []int32) bool { + seen := map[int32]bool{} + for i, a := range c { + if seen[a] { + return false + } + seen[a] = true + for _, b := range c[i+1:] { + if !adj[a][b] { + return false + } + } + } + return true +} + +func TestFindCliqueAgainstBruteForce(t *testing.T) { + rng := rand.New(rand.NewSource(1)) + for trial := 0; trial < 40; trial++ { + n := 8 + rng.Intn(10) + adj := randomGraph(rng, n, 0.3+0.5*rng.Float64()) + maxSize := bruteMaxClique(adj) + filter := matrixFilter(adj) + + for _, workers := range []int{1, 3} { + r := FindClique(context.Background(), n, filter, maxSize, CliqueOptions{Workers: workers}) + if !r.Reached || len(r.Clique) < maxSize || !isClique(adj, r.Clique) { + t.Fatalf("trial %d: target %d not reached: %+v", trial, maxSize, r) + } + + r = FindClique(context.Background(), n, filter, maxSize+1, CliqueOptions{Workers: workers}) + if r.Reached || !r.Exhausted { + t.Fatalf("trial %d: size %d cannot exist, got %+v", trial, maxSize+1, r) + } + // best effort: a valid clique, not necessarily a maximum one + if len(r.Clique) < 1 || len(r.Clique) > maxSize || !isClique(adj, r.Clique) { + t.Fatalf("trial %d: invalid best clique %v (maximum is %d)", trial, r.Clique, maxSize) + } + } + } +} + +func TestFindCliqueTrivialTargets(t *testing.T) { + adj := randomGraph(rand.New(rand.NewSource(2)), 5, 0.5) + if r := FindClique(context.Background(), 5, matrixFilter(adj), 0, CliqueOptions{}); !r.Reached { + t.Fatal("target 0 is always reached") + } + if r := FindClique(context.Background(), 0, matrixFilter(adj), 1, CliqueOptions{}); r.Reached { + t.Fatal("no clique in an empty graph") + } + if r := FindClique(context.Background(), 5, matrixFilter(adj), 1, CliqueOptions{}); !r.Reached || len(r.Clique) != 1 { + t.Fatalf("a single vertex is a clique: %+v", r) + } +} + +// A hopeless search must stop with the deadline and return a valid clique. +// (When the target is far above the clique number of a sparse graph, pruning +// proves the absence at once: see TestFindCliqueAgainstBruteForce.) +func TestFindCliqueTimeout(t *testing.T) { + // dense graph: its largest clique is far below the target, but pruning only + // bites after many levels, so the search cannot finish within the deadline + const n = 300 + adj := randomGraph(rand.New(rand.NewSource(3)), n, 0.9) + ctx, cancel := context.WithTimeout(context.Background(), 300*time.Millisecond) + defer cancel() + + start := time.Now() + r := FindClique(ctx, n, matrixFilter(adj), 150, CliqueOptions{Workers: 2}) + if elapsed := time.Since(start); elapsed > 3*time.Second { + t.Fatalf("deadline not honoured: %v", elapsed) + } + if r.Reached || r.Exhausted { + t.Fatalf("size 150 cannot be reached nor proven absent in time: %+v", r) + } + if len(r.Clique) < 2 || !isClique(adj, r.Clique) { + t.Fatalf("invalid best clique %v", r.Clique) + } +} diff --git a/pkg/obikmer/enumkmer.go b/pkg/obikmer/enumkmer.go index 678bf08..f0c18f5 100644 --- a/pkg/obikmer/enumkmer.go +++ b/pkg/obikmer/enumkmer.go @@ -8,18 +8,32 @@ import "iter" // a negative value means unlimited, 0 forbids the base altogether // (so MaxCount[2] = 0 restricts the alphabet to A, C, T). // +// MaxGC is the maximum number of C or G bases; a negative value means +// unlimited. +// // MaxRun is the maximum number of consecutive identical bases; -// a value <= 0 means unlimited. +// a value <= 0 means unlimited. MinRun requires at least one run of that +// length; a value <= 1 means no requirement. +// +// Accepted and Rejected are optional per-position sets of bases (4 bits masks, +// bit b set when the base of code b is in the set, see ParseIUPACMasks). +// When set, their length must be k. A k-mer is accepted only if each of its +// bases belongs to the Accepted set of its position, and is refused if every +// one of its bases belongs to the Rejected set of its position. // // The zero value forbids every base: use NoKmerConstraints as a starting point. type KmerConstraints struct { MaxCount [4]int + MaxGC int MaxRun int + MinRun int + Accepted []uint8 + Rejected []uint8 } // NoKmerConstraints returns constraints accepting every k-mer. func NoKmerConstraints() KmerConstraints { - return KmerConstraints{MaxCount: [4]int{-1, -1, -1, -1}} + return KmerConstraints{MaxCount: [4]int{-1, -1, -1, -1}, MaxGC: -1} } // IterConstrainedKmers returns an iterator over all k-mers of size k @@ -27,11 +41,14 @@ func NoKmerConstraints() KmerConstraints { // numerical order. // // K-mers are built base by base and a branch is abandoned as soon as a -// constraint is violated (both MaxCount and MaxRun are monotone: a violating -// prefix cannot be repaired by extension), so the cost is proportional to the -// number of explored valid prefixes, not to 4^k. +// constraint is violated or cannot be met any more (MaxCount, MaxGC, MaxRun +// and Accepted are monotone: a violating prefix cannot be repaired by +// extension; MinRun is checked against the positions left), so the cost is +// proportional to the number of explored prefixes, not to 4^k. Rejected is +// checked on complete k-mers. // -// k must be between 1 and 31, otherwise the iterator is empty. +// k must be between 1 and 31, otherwise the iterator is empty. It is also +// empty if Accepted or Rejected is set with a length different from k. // // Example: // @@ -41,7 +58,7 @@ func NoKmerConstraints() KmerConstraints { // for kmer := range IterConstrainedKmers(8, c) { ... } func IterConstrainedKmers(k int, c KmerConstraints) iter.Seq[uint64] { return func(yield func(uint64) bool) { - if k < 1 || k > 31 { + if k < 1 || k > 31 || !c.validPositions(k) { return } @@ -59,16 +76,27 @@ func IterConstrainedKmers(k int, c KmerConstraints) iter.Seq[uint64] { } var counts [4]int + gc := 0 - var walk func(depth int, kmer uint64, last, run int) bool - walk = func(depth int, kmer uint64, last, run int) bool { + var walk func(depth int, kmer uint64, last, run, longest int) bool + walk = func(depth int, kmer uint64, last, run, longest int) bool { if depth == k { + if c.Rejected != nil && matchesAllMasks(kmer, k, c.Rejected) { + return true + } return yield(kmer) } for b := range 4 { if c.MaxCount[b] >= 0 && counts[b] >= c.MaxCount[b] { continue } + if c.Accepted != nil && c.Accepted[depth]>>uint(b)&1 == 0 { + continue + } + isGC := b == 1 || b == 2 + if isGC && c.MaxGC >= 0 && gc >= c.MaxGC { + continue + } r := 1 if b == last { r = run + 1 @@ -76,9 +104,19 @@ func IterConstrainedKmers(k int, c KmerConstraints) iter.Seq[uint64] { if c.MaxRun > 0 && r > c.MaxRun { continue } + l := max(longest, r) + if c.MinRun > 1 && l < c.MinRun && r+(k-depth-1) < c.MinRun { + continue + } counts[b]++ - ok := walk(depth+1, kmer<<2|uint64(b), b, r) + if isGC { + gc++ + } + ok := walk(depth+1, kmer<<2|uint64(b), b, r, l) counts[b]-- + if isGC { + gc-- + } if !ok { return false } @@ -86,6 +124,6 @@ func IterConstrainedKmers(k int, c KmerConstraints) iter.Seq[uint64] { return true } - walk(0, 0, -1, 0) + walk(0, 0, -1, 0, 0) } } diff --git a/pkg/obikmer/kmer_constraints.go b/pkg/obikmer/kmer_constraints.go new file mode 100644 index 0000000..6fcf25c --- /dev/null +++ b/pkg/obikmer/kmer_constraints.go @@ -0,0 +1,141 @@ +package obikmer + +import ( + "fmt" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" +) + +// ParseIUPACMasks converts an IUPAC nucleotide pattern into one 4 bits mask per +// position (bit 0 is A, bit 1 is C, bit 2 is G, bit 3 is T), the format used by +// KmerConstraints.Accepted and Rejected. +func ParseIUPACMasks(pattern string) ([]uint8, error) { + masks := make([]uint8, len(pattern)) + for i := 0; i < len(pattern); i++ { + m := obiseq.IUPACNucMask(pattern[i]) + if m == 0 { + return nil, fmt.Errorf("invalid IUPAC nucleotide code %q at position %d of pattern %q", + pattern[i], i+1, pattern) + } + masks[i] = m + } + return masks, nil +} + +// RestrictAlphabet forbids every base that is not denoted by one of the IUPAC +// codes of alphabet. For instance "ACT" (or "H") forbids G. +func (c *KmerConstraints) RestrictAlphabet(alphabet string) error { + masks, err := ParseIUPACMasks(alphabet) + if err != nil { + return err + } + var allowed uint8 + for _, m := range masks { + allowed |= m + } + if allowed == 0 { + return fmt.Errorf("empty alphabet") + } + for b := range 4 { + if allowed>>uint(b)&1 == 0 { + c.MaxCount[b] = 0 + } + } + return nil +} + +// validPositions tells whether Accepted and Rejected, when set, have length k. +func (c *KmerConstraints) validPositions(k int) bool { + return (c.Accepted == nil || len(c.Accepted) == k) && + (c.Rejected == nil || len(c.Rejected) == k) +} + +// matchesAllMasks tells whether every base of the k-mer belongs to the set +// given for its position. +func matchesAllMasks(kmer uint64, k int, masks []uint8) bool { + for i := k - 1; i >= 0; i-- { + if masks[i]>>(kmer&3)&1 == 0 { + return false + } + kmer >>= 2 + } + return true +} + +// Accepts tells whether a k-mer of size k satisfies every constraint. +// It is the predicate whose solutions IterConstrainedKmers enumerates, to be +// used to filter k-mers coming from elsewhere. +func (c *KmerConstraints) Accepts(kmer uint64, k int) bool { + if k < 1 || k > 31 || !c.validPositions(k) { + return false + } + + var counts [4]int + gc, run, longest, last := 0, 0, 0, -1 + for i := k - 1; i >= 0; i-- { + b := int(kmer>>(2*uint(i))) & 3 + if c.Accepted != nil && c.Accepted[k-1-i]>>uint(b)&1 == 0 { + return false + } + counts[b]++ + if c.MaxCount[b] >= 0 && counts[b] > c.MaxCount[b] { + return false + } + if b == 1 || b == 2 { + gc++ + if c.MaxGC >= 0 && gc > c.MaxGC { + return false + } + } + if b == last { + run++ + } else { + run = 1 + } + last = b + longest = max(longest, run) + if c.MaxRun > 0 && run > c.MaxRun { + return false + } + } + + if c.MinRun > 1 && longest < c.MinRun { + return false + } + if c.Rejected != nil && matchesAllMasks(kmer, k, c.Rejected) { + return false + } + return true +} + +// HammingBound returns the sphere-packing (Hamming) upper bound on the number +// of words of length k over an alphabet of alphabetSize letters that can be +// pairwise at Hamming distance >= d. No larger family exists; a family of +// that size may not exist either. +func HammingBound(alphabetSize, k, d int) float64 { + if d <= 1 { + total := 1.0 + for i := 0; i < k; i++ { + total *= float64(alphabetSize) + } + return total + } + if d%2 == 0 { // A(k, d) <= A(k-1, d-1) + k-- + d-- + } + if k < 1 { + return 1 + } + t := (d - 1) / 2 + total, ball, binom, pow := 1.0, 1.0, 1.0, 1.0 + for i := 0; i < k; i++ { + total *= float64(alphabetSize) + } + for i := 1; i <= t && i <= k; i++ { + binom = binom * float64(k-i+1) / float64(i) + pow *= float64(alphabetSize - 1) + ball += binom * pow + } + return total / ball +} diff --git a/pkg/obikmer/kmer_constraints_test.go b/pkg/obikmer/kmer_constraints_test.go new file mode 100644 index 0000000..0d1f06a --- /dev/null +++ b/pkg/obikmer/kmer_constraints_test.go @@ -0,0 +1,123 @@ +package obikmer + +import ( + "math" + "math/rand" + "testing" +) + +// IterConstrainedKmers must enumerate exactly the k-mers accepted by Accepts. +func TestIterConstrainedKmersMatchesAccepts(t *testing.T) { + rng := rand.New(rand.NewSource(10)) + patterns := []string{"", "nnnnnn", "yrnnnn", "cnnnnn", "annnnt"} + + for trial := 0; trial < 60; trial++ { + k := 6 + c := NoKmerConstraints() + for b := range 4 { + if rng.Intn(3) == 0 { + c.MaxCount[b] = rng.Intn(4) + } + } + if rng.Intn(2) == 0 { + c.MaxGC = rng.Intn(5) + } + if rng.Intn(2) == 0 { + c.MaxRun = 1 + rng.Intn(3) + } + if rng.Intn(2) == 0 { + c.MinRun = 2 + rng.Intn(3) + } + if p := patterns[rng.Intn(len(patterns))]; p != "" { + c.Accepted, _ = ParseIUPACMasks(p) + } + if p := patterns[rng.Intn(len(patterns))]; p != "" && rng.Intn(2) == 0 { + c.Rejected, _ = ParseIUPACMasks(p) + } + + var want []uint64 + for km := uint64(0); km < 1<<(2*uint(k)); km++ { + if c.Accepts(km, k) { + want = append(want, km) + } + } + got := collectConstrained(k, c) + if len(got) != len(want) { + t.Fatalf("trial %d %+v: got %d k-mers, want %d", trial, c, len(got), len(want)) + } + for i := range got { + if got[i] != want[i] { + t.Fatalf("trial %d %+v: mismatch at %d: %d != %d", trial, c, i, got[i], want[i]) + } + } + } +} + +func TestConstraintsGCAndRuns(t *testing.T) { + enc := func(s string) uint64 { return EncodeKmer([]byte(s), len(s)) } + + c := NoKmerConstraints() + c.MaxGC = 2 + if !c.Accepts(enc("ACGTAA"), 6) || c.Accepts(enc("ACGCAA"), 6) { + t.Fatal("MaxGC not enforced") + } + + c = NoKmerConstraints() + c.MinRun = 3 + if !c.Accepts(enc("ACCCGT"), 6) || c.Accepts(enc("ACCGGT"), 6) { + t.Fatal("MinRun not enforced") + } + + c = NoKmerConstraints() + c.Rejected, _ = ParseIUPACMasks("cnnnnn") // no C in first position + if c.Accepts(enc("CAAAAA"), 6) || !c.Accepts(enc("AAAAAA"), 6) { + t.Fatal("Rejected not enforced") + } +} + +func TestRestrictAlphabetAndPatterns(t *testing.T) { + c := NoKmerConstraints() + if err := c.RestrictAlphabet("ACT"); err != nil { + t.Fatal(err) + } + if got := len(collectConstrained(4, c)); got != 81 { + t.Fatalf("expected 3^4 = 81 k-mers, got %d", got) + } + c = NoKmerConstraints() + if err := c.RestrictAlphabet("H"); err != nil { // H = A, C or T + t.Fatal(err) + } + if c.MaxCount != [4]int{-1, -1, 0, -1} { + t.Fatalf("unexpected counts %v", c.MaxCount) + } + if err := c.RestrictAlphabet("AX"); err == nil { + t.Fatal("invalid IUPAC code must be refused") + } + if _, err := ParseIUPACMasks("acgz"); err == nil { + t.Fatal("invalid IUPAC code must be refused") + } + // wrong pattern length + c = NoKmerConstraints() + c.Accepted, _ = ParseIUPACMasks("nnn") + if got := collectConstrained(4, c); len(got) != 0 { + t.Fatalf("pattern of wrong length must give nothing, got %d", len(got)) + } +} + +func TestHammingBound(t *testing.T) { + cases := []struct { + q, k, d int + want float64 + }{ + {4, 7, 3, 16384.0 / 22}, + {4, 8, 4, 16384.0 / 22}, // d even: reduces to (7, 3) + {4, 9, 5, 262144.0 / 352}, + {3, 9, 3, 19683.0 / 19}, + {4, 5, 1, 1024}, + } + for _, c := range cases { + if got := HammingBound(c.q, c.k, c.d); math.Abs(got-c.want) > 1e-6 { + t.Errorf("HammingBound(%d,%d,%d) = %v, want %v", c.q, c.k, c.d, got, c.want) + } + } +} diff --git a/pkg/obikmer/kmer_dissimilarity.go b/pkg/obikmer/kmer_dissimilarity.go new file mode 100644 index 0000000..5f1c57a --- /dev/null +++ b/pkg/obikmer/kmer_dissimilarity.go @@ -0,0 +1,84 @@ +package obikmer + +import ( + "fmt" + "strings" + "sync" +) + +// KmerMetric selects the distance used to compare two k-mers. +type KmerMetric int + +const ( + // HammingMetric counts the positions at which two k-mers differ. + HammingMetric KmerMetric = iota + // LCSMetric is the alignment length minus the LCS length (see KmerLCS). + LCSMetric +) + +// ParseKmerMetric converts "hamming" or "lcs" (case insensitive) to a KmerMetric. +func ParseKmerMetric(name string) (KmerMetric, error) { + switch strings.ToLower(name) { + case "hamming": + return HammingMetric, nil + case "lcs": + return LCSMetric, nil + } + return HammingMetric, fmt.Errorf("unknown distance type %q (expected hamming or lcs)", name) +} + +func (m KmerMetric) String() string { + if m == LCSMetric { + return "lcs" + } + return "hamming" +} + +type lcsScratch struct { + partners []uint64 + results []KmerLCSResult +} + +// NewKmerDissimilarityFilter describes the graph whose vertices are kmers[i] +// (k-mers of size k, 1 <= k <= 31) and where two vertices are linked when +// their distance is at least dmin. The returned function has the signature of +// obigraph.NeighborFilter: it appends to out[:0] the candidates (indices in +// kmers) linked to v. It is safe for concurrent use. +func NewKmerDissimilarityFilter(kmers []uint64, k int, metric KmerMetric, dmin int) func(v int32, candidates, out []int32) []int32 { + if metric == HammingMetric { + return func(v int32, candidates, out []int32) []int32 { + q := kmers[v] + out = out[:0] + for _, c := range candidates { + if HammingDistance(q, kmers[c], k) >= dmin { + out = append(out, c) + } + } + return out + } + } + + pool := sync.Pool{New: func() any { return new(lcsScratch) }} + return func(v int32, candidates, out []int32) []int32 { + s := pool.Get().(*lcsScratch) + defer pool.Put(s) + + s.partners = s.partners[:0] + for _, c := range candidates { + s.partners = append(s.partners, kmers[c]) + } + if cap(s.results) < len(candidates) { + s.results = make([]KmerLCSResult, len(candidates)) + } + s.results = s.results[:len(candidates)] + KmerLCSBatch(kmers[v], s.partners, k, s.results) + + out = out[:0] + for i, c := range candidates { + if s.results[i].Distance() >= dmin { + out = append(out, c) + } + } + return out + } +} diff --git a/pkg/obikmer/kmer_dissimilarity_test.go b/pkg/obikmer/kmer_dissimilarity_test.go new file mode 100644 index 0000000..94a92e6 --- /dev/null +++ b/pkg/obikmer/kmer_dissimilarity_test.go @@ -0,0 +1,113 @@ +package obikmer_test + +import ( + "context" + "math/rand" + "testing" + "time" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obigraph" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer" +) + +func allKmers(k int) []uint64 { + kmers := make([]uint64, 0, 1<<(2*uint(k))) + for km := range obikmer.IterConstrainedKmers(k, obikmer.NoKmerConstraints()) { + kmers = append(kmers, km) + } + return kmers +} + +func TestDissimilarityFilterMatchesDistances(t *testing.T) { + const k = 7 + rng := rand.New(rand.NewSource(1)) + kmers := make([]uint64, 300) + for i := range kmers { + kmers[i] = rng.Uint64() & (1<<(2*k) - 1) + } + candidates := make([]int32, 0, len(kmers)) + for i := 1; i < len(kmers); i++ { + candidates = append(candidates, int32(i)) + } + + for _, metric := range []obikmer.KmerMetric{obikmer.HammingMetric, obikmer.LCSMetric} { + for _, dmin := range []int{1, 3, 5} { + filter := obikmer.NewKmerDissimilarityFilter(kmers, k, metric, dmin) + got := filter(0, candidates, nil) + + var want []int32 + for _, c := range candidates { + var d int + if metric == obikmer.HammingMetric { + d = obikmer.HammingDistance(kmers[0], kmers[c], k) + } else { + d = obikmer.KmerLCS(kmers[0], kmers[c], k).Distance() + } + if d >= dmin { + want = append(want, c) + } + } + if len(got) != len(want) { + t.Fatalf("%v d>=%d: got %d vertices, want %d", metric, dmin, len(got), len(want)) + } + for i := range got { + if got[i] != want[i] { + t.Fatalf("%v d>=%d: mismatch at %d", metric, dmin, i) + } + } + } + } +} + +// A family of 5-mers pairwise at Hamming distance >= 3, built with FindClique. +func TestFindClique_HammingFamily(t *testing.T) { + const k, dmin, target = 5, 3, 32 + kmers := allKmers(k) + filter := obikmer.NewKmerDissimilarityFilter(kmers, k, obikmer.HammingMetric, dmin) + + ctx, cancel := context.WithTimeout(context.Background(), 30*time.Second) + defer cancel() + r := obigraph.FindClique(ctx, len(kmers), filter, target, obigraph.CliqueOptions{Workers: 2}) + if !r.Reached { + t.Fatalf("family of %d not found (best %d, exhausted %v)", target, len(r.Clique), r.Exhausted) + } + for i, a := range r.Clique { + for _, b := range r.Clique[i+1:] { + if obikmer.HammingDistance(kmers[a], kmers[b], k) < dmin { + t.Fatalf("%d and %d are too close", kmers[a], kmers[b]) + } + } + } +} + +func TestFindClique_LCSFamily(t *testing.T) { + const k, dmin, target = 6, 3, 20 + kmers := allKmers(k) + filter := obikmer.NewKmerDissimilarityFilter(kmers, k, obikmer.LCSMetric, dmin) + + ctx, cancel := context.WithTimeout(context.Background(), 60*time.Second) + defer cancel() + r := obigraph.FindClique(ctx, len(kmers), filter, target, obigraph.CliqueOptions{Workers: 2}) + if !r.Reached { + t.Fatalf("family of %d not found (best %d, exhausted %v)", target, len(r.Clique), r.Exhausted) + } + for i, a := range r.Clique { + for _, b := range r.Clique[i+1:] { + if obikmer.KmerLCS(kmers[a], kmers[b], k).Distance() < dmin { + t.Fatalf("%d and %d are too close", kmers[a], kmers[b]) + } + } + } +} + +func TestParseKmerMetric(t *testing.T) { + if m, err := obikmer.ParseKmerMetric("LCS"); err != nil || m != obikmer.LCSMetric { + t.Fatal("lcs not parsed") + } + if m, err := obikmer.ParseKmerMetric("hamming"); err != nil || m != obikmer.HammingMetric { + t.Fatal("hamming not parsed") + } + if _, err := obikmer.ParseKmerMetric("levenshtein"); err == nil { + t.Fatal("unknown metric must be refused") + } +} diff --git a/pkg/obiseq/iupac_nog.go b/pkg/obiseq/iupac_nog.go index a3fba9b..7565263 100644 --- a/pkg/obiseq/iupac_nog.go +++ b/pkg/obiseq/iupac_nog.go @@ -26,3 +26,16 @@ func SameIUPACNuc(a, b byte) bool { } return a == b } + +// IUPACNucMask returns the set of nucleotides denoted by the IUPAC code b +// (case insensitive, U is read as T) as a 4 bits mask: bit 0 is A, bit 1 is C, +// bit 2 is G and bit 3 is T. It returns 0 if b is not an IUPAC nucleotide code. +func IUPACNucMask(b byte) byte { + if (b >= 'A') && (b <= 'Z') { + b |= 32 + } + if (b < 'a') || (b > 'z') { + return 0 + } + return _iupac[b-'a'] +} diff --git a/pkg/obitools/obioligotag/input.go b/pkg/obitools/obioligotag/input.go new file mode 100644 index 0000000..992ba8e --- /dev/null +++ b/pkg/obitools/obioligotag/input.go @@ -0,0 +1,148 @@ +package obioligotag + +import ( + "bufio" + "fmt" + "io" + "slices" + "strings" + + log "github.com/sirupsen/logrus" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obiconvert" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiutils" +) + +// maxOligoSize is the longest oligonucleotide handled. +const maxOligoSize = 12 + +// isPlainOligoLine tells whether a line looks like an oligonucleotide: only +// nucleotide letters, which a sequence file title line or a header never is. +func isPlainOligoLine(line string) bool { + if line == "" { + return false + } + for i := 0; i < len(line); i++ { + if obiseq.IUPACNucMask(line[i]) == 0 { + return false + } + } + return true +} + +// isPlainOligoFile guesses whether filename is a text file with one +// oligonucleotide per line, by looking at its first significant line. +func isPlainOligoFile(filename string) (bool, error) { + file, err := obiutils.Ropen(filename) + if err != nil { + return false, err + } + defer file.Close() + + reader := bufio.NewReader(file) + for { + line, err := reader.ReadString('\n') + line = strings.TrimSpace(line) + if line != "" && !strings.HasPrefix(line, "#") { + return isPlainOligoLine(line), nil + } + if err != nil { + if err == io.EOF { + return false, fmt.Errorf("file %s is empty", filename) + } + return false, err + } + } +} + +// oligoCollector accumulates encoded oligonucleotides of a single size. +type oligoCollector struct { + size int + kmers []uint64 + skipped int +} + +func (oc *oligoCollector) add(seq []byte) error { + if len(seq) == 0 { + return nil + } + if len(seq) > maxOligoSize { + return fmt.Errorf("oligonucleotide %q is longer than the maximum size %d", seq, maxOligoSize) + } + if oc.size == 0 { + oc.size = len(seq) + } + if len(seq) != oc.size { + return fmt.Errorf("all the oligonucleotides must have the same length: found %d and %d", oc.size, len(seq)) + } + for _, b := range seq { + switch b | 32 { + case 'a', 'c', 'g', 't', 'u': + default: + oc.skipped++ + return nil + } + } + oc.kmers = append(oc.kmers, obikmer.EncodeKmer(seq, len(seq))) + return nil +} + +// readOligoList reads the oligonucleotides of filename, either as a sequence +// file or as a text file with one oligonucleotide per line. Oligonucleotides +// containing something else than A, C, G, T or U are skipped, and +// duplicates removed. The result is sorted. +func readOligoList(filename string) (kmers []uint64, size int, err error) { + plain := false + if obiconvert.CLIInputFormat() == "guessed" { + plain, err = isPlainOligoFile(filename) + if err != nil { + return nil, 0, err + } + } + + oc := &oligoCollector{} + + if plain { + file, err := obiutils.Ropen(filename) + if err != nil { + return nil, 0, err + } + defer file.Close() + + scanner := bufio.NewScanner(file) + for scanner.Scan() { + line := strings.TrimSpace(scanner.Text()) + if line == "" || strings.HasPrefix(line, "#") { + continue + } + if err := oc.add([]byte(line)); err != nil { + return nil, 0, err + } + } + if err := scanner.Err(); err != nil { + return nil, 0, err + } + } else { + sequences, err := obiconvert.CLIReadBioSequences(filename) + if err != nil { + return nil, 0, err + } + for sequences.Next() { + for _, s := range sequences.Get().Slice() { + if err := oc.add(s.Sequence()); err != nil { + return nil, 0, err + } + } + } + } + + if oc.skipped > 0 { + log.Warnf("%d oligonucleotides containing characters other than A, C, G, T were ignored", oc.skipped) + } + + slices.Sort(oc.kmers) + oc.kmers = slices.Compact(oc.kmers) + return oc.kmers, oc.size, nil +} diff --git a/pkg/obitools/obioligotag/obioligotag.go b/pkg/obitools/obioligotag/obioligotag.go new file mode 100644 index 0000000..bb3e93b --- /dev/null +++ b/pkg/obitools/obioligotag/obioligotag.go @@ -0,0 +1,217 @@ +package obioligotag + +import ( + "context" + "fmt" + "math" + "slices" + "time" + + log "github.com/sirupsen/logrus" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obigraph" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiiter" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" +) + +// buildConstraints translates the command line options into constraints on +// oligonucleotides of the given size. +func buildConstraints(size int) (obikmer.KmerConstraints, error) { + c := obikmer.NoKmerConstraints() + + if err := c.RestrictAlphabet(CLIAlphabet()); err != nil { + return c, fmt.Errorf("invalid --alphabet: %w", err) + } + + c.MaxGC = CLIGCMax() + if CLIHomopolymerMax() > 0 { + c.MaxRun = CLIHomopolymerMax() + } + c.MinRun = CLIHomopolymerMin() + + if pattern := CLIAccepted(); pattern != "" { + if len(pattern) != size { + return c, fmt.Errorf("--accepted pattern length (%d) differs from the oligonucleotide size (%d)", len(pattern), size) + } + masks, err := obikmer.ParseIUPACMasks(pattern) + if err != nil { + return c, fmt.Errorf("invalid --accepted: %w", err) + } + c.Accepted = masks + } + + if pattern := CLIRejected(); pattern != "" { + if len(pattern) != size { + return c, fmt.Errorf("--rejected pattern length (%d) differs from the oligonucleotide size (%d)", len(pattern), size) + } + masks, err := obikmer.ParseIUPACMasks(pattern) + if err != nil { + return c, fmt.Errorf("invalid --rejected: %w", err) + } + c.Rejected = masks + } + + return c, nil +} + +// candidateOligos returns the sorted candidate oligonucleotides and their size: +// either every oligonucleotide satisfying the constraints, or those of the +// --oligo-list file satisfying them. +func candidateOligos() ([]uint64, int, error) { + size := CLIOligoSize() + list := CLIOligoList() + + switch { + case list != "" && size != 0: + return nil, 0, fmt.Errorf("--oligo-size and --oligo-list cannot be used together") + case list == "" && size == 0: + return nil, 0, fmt.Errorf("either --oligo-size or --oligo-list is required") + case list == "" && (size < 1 || size > maxOligoSize): + return nil, 0, fmt.Errorf("--oligo-size must be between 1 and %d", maxOligoSize) + } + + var kmers []uint64 + if list != "" { + var err error + kmers, size, err = readOligoList(list) + if err != nil { + return nil, 0, fmt.Errorf("cannot read %s: %w", list, err) + } + log.Infof("%d oligonucleotides of size %d read from %s", len(kmers), size, list) + } + + constraints, err := buildConstraints(size) + if err != nil { + return nil, 0, err + } + + if list == "" { + kmers = slices.Collect(obikmer.IterConstrainedKmers(size, constraints)) + } else { + kmers = slices.DeleteFunc(kmers, func(km uint64) bool { return !constraints.Accepts(km, size) }) + } + + if len(kmers) == 0 { + return nil, 0, fmt.Errorf("no oligonucleotide satisfies the constraints") + } + log.Infof("%d candidate oligonucleotides of size %d", len(kmers), size) + return kmers, size, nil +} + +// checkFeasibility stops the program early when the family size is beyond what +// the Hamming bound allows. +func checkFeasibility(size int, metric obikmer.KmerMetric) error { + if metric != obikmer.HammingMetric { + return nil + } + alphabet := obikmer.NoKmerConstraints() + if err := alphabet.RestrictAlphabet(CLIAlphabet()); err != nil { + return err + } + q := 0 + for _, m := range alphabet.MaxCount { + if m != 0 { + q++ + } + } + if CLIOligoList() != "" { + q = 4 // the candidates of a list are not restricted by the alphabet of the generator + } + bound := obikmer.HammingBound(q, size, CLIDistance()) + if float64(CLIFamilySize()) > math.Floor(bound+1e-9) { + return fmt.Errorf("no family of %d oligonucleotides of size %d at Hamming distance >= %d exists "+ + "(sphere-packing bound: %.0f)", CLIFamilySize(), size, CLIDistance(), math.Floor(bound+1e-9)) + } + return nil +} + +// CLIOligoFamily builds the family of oligonucleotides described by the +// command line options: it returns oligonucleotides, sorted, at a distance of at +// least --distance from each other. If no family of the required size is +// found before the time limit, the largest family found is returned together +// with a warning. +func CLIOligoFamily() ([]string, error) { + target := CLIFamilySize() + dmin := CLIDistance() + + if target < 1 { + return nil, fmt.Errorf("--family-size is required and must be at least 1") + } + if dmin < 1 { + return nil, fmt.Errorf("--distance must be at least 1") + } + metric, err := obikmer.ParseKmerMetric(CLIDistanceType()) + if err != nil { + return nil, err + } + + kmers, size, err := candidateOligos() + if err != nil { + return nil, err + } + if err := checkFeasibility(size, metric); err != nil { + return nil, err + } + + ctx := context.Background() + if CLITimeout() > 0 { + var cancel context.CancelFunc + ctx, cancel = context.WithTimeout(ctx, time.Duration(CLITimeout())*time.Second) + defer cancel() + } + + filter := obikmer.NewKmerDissimilarityFilter(kmers, size, metric, dmin) + result := obigraph.FindClique(ctx, len(kmers), filter, target, + obigraph.CliqueOptions{Workers: obidefault.ParallelWorkers(), Seed: 1}) + + switch { + case result.Reached: + log.Infof("family of %d oligonucleotides found", len(result.Clique)) + case result.Exhausted: + log.Warnf("no family of %d oligonucleotides exists with these constraints: "+ + "the best family found has %d oligonucleotides", target, len(result.Clique)) + default: + log.Warnf("time limit reached before finding a family of %d oligonucleotides: "+ + "the best family found has %d oligonucleotides", target, len(result.Clique)) + } + + family := make([]uint64, len(result.Clique)) + for i, v := range result.Clique { + family[i] = kmers[v] + } + slices.Sort(family) + + oligos := make([]string, len(family)) + for i, km := range family { + oligos[i] = string(obikmer.DecodeKmer(km, size, nil)) + } + return oligos, nil +} + +// CLIOligosToIterator wraps oligonucleotides in a sequence iterator, with the +// identifiers oligo_001, oligo_002... +func CLIOligosToIterator(oligos []string) obiiter.IBioSequence { + width := len(fmt.Sprint(len(oligos))) + batchSize := max(obidefault.BatchSize(), 1) + + iterator := obiiter.MakeIBioSequence() + iterator.Add(1) + + go func() { + defer iterator.Done() + for order, start := 0, 0; start < len(oligos); order, start = order+1, start+batchSize { + end := min(start+batchSize, len(oligos)) + slice := make(obiseq.BioSequenceSlice, 0, end-start) + for i := start; i < end; i++ { + slice = append(slice, + obiseq.NewBioSequence(fmt.Sprintf("oligo_%0*d", width, i+1), []byte(oligos[i]), "")) + } + iterator.Push(obiiter.MakeBioSequenceBatch("obioligotag", order, slice)) + } + }() + + go iterator.WaitAndClose() + return iterator +} diff --git a/pkg/obitools/obioligotag/obioligotag_test.go b/pkg/obitools/obioligotag/obioligotag_test.go new file mode 100644 index 0000000..5310d07 --- /dev/null +++ b/pkg/obitools/obioligotag/obioligotag_test.go @@ -0,0 +1,177 @@ +package obioligotag + +import ( + "os" + "path/filepath" + "strings" + "testing" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer" +) + +// defaults puts every option back to its default value. +func defaults() { + SetOligoSize(0) + SetOligoList("") + SetFamilySize(0) + SetDistance(1) + SetDistanceType("hamming") + SetAlphabet("ACGT") + SetGCMax(-1) + SetAccepted("") + SetRejected("") + SetHomopolymerMax(-1) + SetHomopolymerMin(0) + SetTimeout(60) +} + +// reset starts a test from the default options and restores them at its end. +func reset(t *testing.T) { + t.Helper() + defaults() + t.Cleanup(defaults) +} + +func hamming(a, b string) int { + d := 0 + for i := range a { + if a[i] != b[i] { + d++ + } + } + return d +} + +func TestOligoFamilyHamming(t *testing.T) { + reset(t) + SetOligoSize(5) + SetFamilySize(24) + SetDistance(3) + SetGCMax(3) + SetHomopolymerMax(2) + + family, err := CLIOligoFamily() + if err != nil { + t.Fatal(err) + } + if len(family) < 24 { + t.Fatalf("expected at least 24 oligonucleotides, got %d", len(family)) + } + for i, a := range family { + if strings.Count(a, "g")+strings.Count(a, "c") > 3 { + t.Errorf("%s has too many G/C", a) + } + for _, h := range []string{"aaa", "ccc", "ggg", "ttt"} { + if strings.Contains(a, h) { + t.Errorf("%s has a homopolymer longer than 2", a) + } + } + for _, b := range family[i+1:] { + if hamming(a, b) < 3 { + t.Errorf("%s and %s are too close", a, b) + } + } + } +} + +func TestOligoFamilyLCSWithoutG(t *testing.T) { + reset(t) + SetOligoSize(6) + SetFamilySize(20) + SetDistance(3) + SetDistanceType("lcs") + SetAlphabet("ACT") + + family, err := CLIOligoFamily() + if err != nil { + t.Fatal(err) + } + if len(family) < 20 { + t.Fatalf("expected at least 20 oligonucleotides, got %d", len(family)) + } + for i, a := range family { + if strings.Contains(a, "g") { + t.Errorf("%s contains a G", a) + } + for _, b := range family[i+1:] { + ka, kb := obikmer.EncodeKmer([]byte(a), 6), obikmer.EncodeKmer([]byte(b), 6) + if obikmer.KmerLCS(ka, kb, 6).Distance() < 3 { + t.Errorf("%s and %s are too close", a, b) + } + } + } +} + +func TestOligoListPlainTextAndFasta(t *testing.T) { + dir := t.TempDir() + plain := filepath.Join(dir, "oligos.txt") + fasta := filepath.Join(dir, "oligos.fasta") + os.WriteFile(plain, []byte("# tags\nacgtac\nacgtat\nttttgg\nggcaac\ncaccat\n\nacgtac\nacgNac\n"), 0o644) + os.WriteFile(fasta, []byte(">a\nACGTAC\n>b\nACGTAT\n>c\nTTTTGG\n>d\nGGCAAC\n>e\nCACCAT\n"), 0o644) + + for _, file := range []string{plain, fasta} { + reset(t) + SetOligoList(file) + SetFamilySize(3) + SetDistance(3) + + kmers, size, err := candidateOligos() + if err != nil { + t.Fatalf("%s: %v", file, err) + } + if size != 6 || len(kmers) != 5 { // duplicates removed, 'acgNac' ignored in the text file + t.Fatalf("%s: got %d oligonucleotides of size %d", file, len(kmers), size) + } + family, err := CLIOligoFamily() + if err != nil || len(family) < 3 { + t.Fatalf("%s: family %v, error %v", file, family, err) + } + } +} + +func TestOligoListFilteredByConstraints(t *testing.T) { + reset(t) + file := filepath.Join(t.TempDir(), "oligos.txt") + os.WriteFile(file, []byte("aaaaaa\nacgtac\nacgtat\n"), 0o644) + SetOligoList(file) + SetFamilySize(1) + SetHomopolymerMax(2) + + kmers, _, err := candidateOligos() + if err != nil || len(kmers) != 2 { + t.Fatalf("expected the homopolymer to be filtered out, got %d (%v)", len(kmers), err) + } +} + +func TestOligoTagErrors(t *testing.T) { + cases := map[string]func(){ + "no size nor list": func() { SetFamilySize(2) }, + "size and list": func() { SetFamilySize(2); SetOligoSize(5); SetOligoList("x") }, + "no family size": func() { SetOligoSize(5) }, + "size too large": func() { SetFamilySize(2); SetOligoSize(13) }, + "bad alphabet": func() { SetFamilySize(2); SetOligoSize(5); SetAlphabet("ACX") }, + "bad pattern length": func() { SetFamilySize(2); SetOligoSize(5); SetAccepted("nnn") }, + "bad distance type": func() { SetFamilySize(2); SetOligoSize(5); SetDistanceType("levenshtein") }, + "beyond Hamming bound": func() { SetFamilySize(100); SetOligoSize(5); SetDistance(3) }, + "empty candidate set": func() { SetFamilySize(2); SetOligoSize(4); SetAlphabet("A"); SetHomopolymerMax(2) }, + } + for name, setup := range cases { + t.Run(name, func(t *testing.T) { + reset(t) + setup() + if family, err := CLIOligoFamily(); err == nil { + t.Errorf("expected an error, got %v", family) + } + }) + } +} + +func TestOligoIsPlainLine(t *testing.T) { + for line, want := range map[string]bool{ + "acgt": true, "ACGTRYN": true, ">seq1": false, "@read": false, "a,b": false, "": false, + } { + if got := isPlainOligoLine(line); got != want { + t.Errorf("isPlainOligoLine(%q) = %v, want %v", line, got, want) + } + } +} diff --git a/pkg/obitools/obioligotag/options.go b/pkg/obitools/obioligotag/options.go new file mode 100644 index 0000000..4ba194a --- /dev/null +++ b/pkg/obitools/obioligotag/options.go @@ -0,0 +1,150 @@ +// obioligotag function utility package. +// +// The obitools/obioligotag package contains every +// functions specificaly required by the obioligotag utility. +package obioligotag + +import ( + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obiconvert" + "github.com/DavidGamba/go-getoptions" +) + +var __oligo_size__ = 0 +var __oligo_list__ = "" +var __family_size__ = 0 +var __distance__ = 1 +var __distance_type__ = "hamming" +var __alphabet__ = "ACGT" +var __gc_max__ = -1 +var __accepted__ = "" +var __rejected__ = "" +var __homopolymer__ = -1 +var __homopolymer_min__ = 0 +var __timeout__ = 60 + +// OligoTagOptionSet registers the options specific to obioligotag. +func OligoTagOptionSet(options *getoptions.GetOpt) { + options.IntVar(&__oligo_size__, "oligo-size", __oligo_size__, + options.Alias("s"), + options.ArgName("SIZE"), + options.Description("Length of the oligonucleotides to generate (1 to 12). "+ + "Cannot be combined with --oligo-list.")) + + options.StringVar(&__oligo_list__, "oligo-list", __oligo_list__, + options.Alias("L"), + options.ArgName("FILENAME"), + options.Description("File of candidate oligonucleotides, in any sequence file format "+ + "or as a text file with one oligonucleotide per line (format detected automatically). "+ + "The family is selected among them. Cannot be combined with --oligo-size.")) + + options.IntVar(&__family_size__, "family-size", __family_size__, + options.Alias("f"), + options.ArgName("SIZE"), + options.Description("Minimum number of oligonucleotides of the family to build.")) + + options.IntVar(&__distance__, "distance", __distance__, + options.Alias("d"), + options.ArgName("DISTANCE"), + options.Description("Minimum distance between any two oligonucleotides of the family.")) + + options.StringVar(&__distance_type__, "distance-type", __distance_type__, + options.ArgName("TYPE"), + options.Description("Distance between oligonucleotides: 'hamming' (number of differing positions) "+ + "or 'lcs' (alignment length minus the length of the longest common subsequence).")) + + options.StringVar(&__alphabet__, "alphabet", __alphabet__, + options.ArgName("LETTERS"), + options.Description("Nucleotides allowed in the oligonucleotides, as IUPAC codes "+ + "(for instance ACT, or H, to forbid G).")) + + options.IntVar(&__gc_max__, "gc-max", __gc_max__, + options.Alias("g"), + options.ArgName("COUNT"), + options.Description("Maximum number of G or C per oligonucleotide.")) + + options.StringVar(&__accepted__, "accepted", __accepted__, + options.Alias("a"), + options.ArgName("IUPAC_PATTERN"), + options.Description("Selected oligonucleotides must match this IUPAC pattern "+ + "(of the same length as the oligonucleotides).")) + + options.StringVar(&__rejected__, "rejected", __rejected__, + options.Alias("r"), + options.ArgName("IUPAC_PATTERN"), + options.Description("Selected oligonucleotides must not match this IUPAC pattern "+ + "(of the same length as the oligonucleotides).")) + + options.IntVar(&__homopolymer__, "homopolymer", __homopolymer__, + options.Alias("p"), + options.ArgName("LENGTH"), + options.Description("Excludes the oligonucleotides containing a homopolymer longer than LENGTH.")) + + options.IntVar(&__homopolymer_min__, "homopolymer-min", __homopolymer_min__, + options.Alias("P"), + options.ArgName("LENGTH"), + options.Description("Keeps only the oligonucleotides containing a homopolymer of at least LENGTH.")) + + options.IntVar(&__timeout__, "timeout", __timeout__, + options.Alias("T"), + options.ArgName("SECONDS"), + options.Description("Time limit, in seconds, for finding a family of the required size "+ + "(0 for no limit). When it is reached, the best family found is returned with a warning.")) +} + +// OptionSet registers every option of obioligotag. +func OptionSet(options *getoptions.GetOpt) { + OligoTagOptionSet(options) + obiconvert.InputOptionSet(options) + obiconvert.OutputOptionSet(options) +} + +// CLIOligoSize returns the size of the oligonucleotides to generate. +func CLIOligoSize() int { return __oligo_size__ } + +// CLIOligoList returns the name of the file of candidate oligonucleotides. +func CLIOligoList() string { return __oligo_list__ } + +// CLIFamilySize returns the minimum size of the family to build. +func CLIFamilySize() int { return __family_size__ } + +// CLIDistance returns the minimum distance between two oligonucleotides. +func CLIDistance() int { return __distance__ } + +// CLIDistanceType returns the name of the distance used. +func CLIDistanceType() string { return __distance_type__ } + +// CLIAlphabet returns the IUPAC codes of the allowed nucleotides. +func CLIAlphabet() string { return __alphabet__ } + +// CLIGCMax returns the maximum number of G or C (negative if unlimited). +func CLIGCMax() int { return __gc_max__ } + +// CLIAccepted returns the IUPAC pattern the oligonucleotides must match. +func CLIAccepted() string { return __accepted__ } + +// CLIRejected returns the IUPAC pattern the oligonucleotides must not match. +func CLIRejected() string { return __rejected__ } + +// CLIHomopolymerMax returns the longest homopolymer allowed (<= 0 if unlimited). +func CLIHomopolymerMax() int { return __homopolymer__ } + +// CLIHomopolymerMin returns the length of the homopolymer required (<= 1 if none). +func CLIHomopolymerMin() int { return __homopolymer_min__ } + +// CLITimeout returns the time limit in seconds (0 if none). +func CLITimeout() int { return __timeout__ } + +// Setters, mainly for testing. + +func SetOligoSize(v int) { __oligo_size__ = v } +func SetOligoList(v string) { __oligo_list__ = v } +func SetFamilySize(v int) { __family_size__ = v } +func SetDistance(v int) { __distance__ = v } +func SetDistanceType(v string) { __distance_type__ = v } +func SetAlphabet(v string) { __alphabet__ = v } +func SetGCMax(v int) { __gc_max__ = v } +func SetAccepted(v string) { __accepted__ = v } +func SetRejected(v string) { __rejected__ = v } +func SetHomopolymerMax(v int) { __homopolymer__ = v } +func SetHomopolymerMin(v int) { __homopolymer_min__ = v } +func SetTimeout(v int) { __timeout__ = v }