Release 4.5.0

Corrects expected values for division and comparison operations in the uint128 test suite. Updates obilandmark to use log.Fatalf instead of log.Fatal, ensuring sequence ID, taxid, and taxonomy name are correctly interpolated in fatal error messages.

.github/workflows/obitest.yml

@@ -10,10 +10,10 @@ jobs:
     runs-on: ubuntu-latest
     steps:
     - name: Setup Go
-      uses: actions/setup-go@v2
+      uses: actions/setup-go@v5
       with:
         go-version: '1.23'
     - name: Checkout obitools4 project
-      uses: actions/checkout@v4
+      uses: actions/checkout@v5
     - name: Run tests
       run: make githubtests

.github/workflows/release.yml

@@ -16,9 +16,9 @@ jobs:
       - name: Setup Go
         uses: actions/setup-go@v5
         with:
-          go-version: "1.23"
+          go-version: "1.26"
       - name: Checkout obitools4 project
-        uses: actions/checkout@v4
+        uses: actions/checkout@v5
       - name: Run tests
         run: make githubtests
 
@@ -49,12 +49,12 @@ jobs:
 
     steps:
       - name: Checkout code
-        uses: actions/checkout@v4
+        uses: actions/checkout@v5
 
       - name: Setup Go
         uses: actions/setup-go@v5
         with:
-          go-version: "1.23"
+          go-version: "1.26"
 
       - name: Extract version from tag
         id: get_version
@@ -62,12 +62,6 @@ jobs:
           TAG=${GITHUB_REF#refs/tags/Release_}
           echo "version=$TAG" >> $GITHUB_OUTPUT
 
-      - name: Install build tools (Linux)
-        if: runner.os == 'Linux'
-        run: |
-          sudo apt-get update -q
-          sudo apt-get install -y musl-tools zlib1g-dev
-
       - name: Install build tools (macOS)
         if: runner.os == 'macOS'
         run: |
@@ -75,21 +69,30 @@ jobs:
           xcode-select --install 2>/dev/null || true
           xcode-select -p
 
-      - name: Build binaries
+      - name: Build binaries (Linux)
+        if: runner.os == 'Linux'
+        env:
+          VERSION: ${{ steps.get_version.outputs.version }}
+        run: |
+          docker run --rm \
+            -v "$(pwd):/src" \
+            -w /src \
+            -e VERSION="${VERSION}" \
+            golang:1.26-alpine \
+            sh -c "apk add --no-cache gcc musl-dev zlib-dev zlib-static make && \
+                   make LDFLAGS='-linkmode=external -extldflags=-static' obitools"
+          mkdir -p artifacts
+          tar -czf artifacts/obitools4_${VERSION}_${{ matrix.output_name }}.tar.gz -C build .
+
+      - name: Build binaries (macOS)
+        if: runner.os == 'macOS'
         env:
           GOOS: ${{ matrix.goos }}
           GOARCH: ${{ matrix.goarch }}
           VERSION: ${{ steps.get_version.outputs.version }}
-          CC: ${{ matrix.goos == 'linux' && 'musl-gcc' || '' }}
-          CGO_CFLAGS: ${{ matrix.goos == 'linux' && '-I/usr/include' || '' }}
         run: |
-          if [ "$GOOS" = "linux" ]; then
-            make LDFLAGS='-linkmode=external -extldflags=-static' obitools
-          else
           make obitools
-          fi
           mkdir -p artifacts
-          # Create a single tar.gz with all binaries for this platform
           tar -czf artifacts/obitools4_${VERSION}_${{ matrix.output_name }}.tar.gz -C build .
 
       - name: Upload artifacts
@@ -104,7 +107,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout code
-        uses: actions/checkout@v4
+        uses: actions/checkout@v5
         with:
           fetch-depth: 0
 

.gitignore

@@ -21,9 +21,11 @@ xx
 
 .rhistory
 /.vscode
+
+/benchmarck
 /build
 /bugs
-
+autodoc
 /ncbitaxo
 
 !/obitests/**

.serena/.gitignore

@@ -0,0 +1,2 @@
+/cache
+/project.local.yml

.serena/project.yml

@@ -0,0 +1,133 @@
+# the name by which the project can be referenced within Serena
+project_name: "obitools4"
+
+
+# list of languages for which language servers are started; choose from:
+#   al                  angular             ansible             bash                clojure
+#   cpp                 cpp_ccls            crystal             csharp              csharp_omnisharp
+#   dart                elixir              elm                 erlang              fortran
+#   fsharp              go                  groovy              haskell             haxe
+#   hlsl                html                java                json                julia
+#   kotlin              lean4               lua                 luau                markdown
+#   matlab              msl                 nix                 ocaml               pascal
+#   perl                php                 php_phpactor        powershell          python
+#   python_jedi         python_ty           r                   rego                ruby
+#   ruby_solargraph     rust                scala               scss                solidity
+#   svelte              swift               systemverilog       terraform           toml
+#   typescript          typescript_vts      vue                 yaml                zig
+#   (This list may be outdated. For the current list, see values of Language enum here:
+#   https://github.com/oraios/serena/blob/main/src/solidlsp/ls_config.py
+#   For some languages, there are alternative language servers, e.g. csharp_omnisharp, ruby_solargraph.)
+# Note:
+#   - For C, use cpp
+#   - For JavaScript, use typescript
+#   - For Angular projects, use angular (subsumes typescript+html; requires `npm install` in the project root)
+#   - For Svelte projects, use svelte (subsumes typescript/javascript for .svelte projects; requires npm)
+#   - For SCSS / Sass / plain CSS, use scss (some-sass-language-server handles all three)
+#   - For Free Pascal/Lazarus, use pascal
+# Special requirements:
+#   Some languages require additional setup/installations.
+#   See here for details: https://oraios.github.io/serena/01-about/020_programming-languages.html#language-servers
+# When using multiple languages, the first language server that supports a given file will be used for that file.
+# The first language is the default language and the respective language server will be used as a fallback.
+# Note that when using the JetBrains backend, language servers are not used and this list is correspondingly ignored.
+languages:
+- go
+
+# the encoding used by text files in the project
+# For a list of possible encodings, see https://docs.python.org/3.11/library/codecs.html#standard-encodings
+encoding: "utf-8"
+
+# line ending convention to use when writing source files.
+# Possible values: unset (use global setting), "lf", "crlf", or "native" (platform default)
+# This does not affect Serena's own files (e.g. memories and configuration files), which always use native line endings.
+line_ending:
+
+# The language backend to use for this project.
+# If not set, the global setting from serena_config.yml is used.
+# Valid values: LSP, JetBrains
+# Note: the backend is fixed at startup. If a project with a different backend
+# is activated post-init, an error will be returned.
+language_backend:
+
+# whether to use project's .gitignore files to ignore files
+ignore_all_files_in_gitignore: true
+
+# advanced configuration option allowing to configure language server-specific options.
+# Maps the language key to the options.
+# Have a look at the docstring of the constructors of the LS implementations within solidlsp (e.g., for C# or PHP) to see which options are available.
+# No documentation on options means no options are available.
+ls_specific_settings: {}
+
+# list of additional workspace folder paths for cross-package reference support (e.g. in monorepos).
+# Paths can be absolute or relative to the project root.
+# Each folder is registered as an LSP workspace folder, enabling language servers to discover
+# symbols and references across package boundaries.
+# Currently supported for: TypeScript.
+# Example:
+#   additional_workspace_folders:
+#     - ../sibling-package
+#     - ../shared-lib
+additional_workspace_folders: []
+
+# list of additional paths to ignore in this project.
+# Same syntax as gitignore, so you can use * and **.
+# Note: global ignored_paths from serena_config.yml are also applied additively.
+ignored_paths: []
+
+# whether the project is in read-only mode
+# If set to true, all editing tools will be disabled and attempts to use them will result in an error
+# Added on 2025-04-18
+read_only: false
+
+# list of tool names to exclude.
+# This extends the existing exclusions (e.g. from the global configuration)
+# Find the list of tools here: https://oraios.github.io/serena/01-about/035_tools.html
+excluded_tools: []
+
+# list of tools to include that would otherwise be disabled (particularly optional tools that are disabled by default).
+# This extends the existing inclusions (e.g. from the global configuration).
+# Find the list of tools here: https://oraios.github.io/serena/01-about/035_tools.html
+included_optional_tools: []
+
+# fixed set of tools to use as the base tool set (if non-empty), replacing Serena's default set of tools.
+# This cannot be combined with non-empty excluded_tools or included_optional_tools.
+# Find the list of tools here: https://oraios.github.io/serena/01-about/035_tools.html
+fixed_tools: []
+
+# list of mode names that are to be activated by default, overriding the setting in the global configuration.
+# The full set of modes to be activated is base_modes (from global config) + default_modes + added_modes.
+# If the setting is undefined/empty, the default_modes from the global configuration (serena_config.yml) apply.
+# Otherwise, this overrides the setting from the global configuration (serena_config.yml).
+# Therefore, you can set this to [] if you do not want the default modes defined in the global config to apply
+# for this project.
+# This setting can, in turn, be overridden by CLI parameters (--mode).
+# See https://oraios.github.io/serena/02-usage/050_configuration.html#modes
+default_modes:
+
+# list of mode names to be activated additionally for this project, e.g. ["query-projects"]
+# The full set of modes to be activated is base_modes (from global config) + default_modes + added_modes.
+# See https://oraios.github.io/serena/02-usage/050_configuration.html#modes
+added_modes:
+
+# initial prompt for the project. It will always be given to the LLM upon activating the project
+# (contrary to the memories, which are loaded on demand).
+initial_prompt: ""
+
+# time budget (seconds) per tool call for the retrieval of additional symbol information
+# such as docstrings or parameter information.
+# This overrides the corresponding setting in the global configuration; see the documentation there.
+# If null or missing, use the setting from the global configuration.
+symbol_info_budget:
+
+# list of regex patterns which, when matched, mark a memory entry as read‑only.
+# Extends the list from the global configuration, merging the two lists.
+read_only_memory_patterns: []
+
+# list of regex patterns for memories to completely ignore.
+# Matching memories will not appear in list_memories or activate_project output
+# and cannot be accessed via read_memory or write_memory.
+# To access ignored memory files, use the read_file tool on the raw file path.
+# Extends the list from the global configuration, merging the two lists.
+# Example: ["_archive/.*", "_episodes/.*"]
+ignored_memory_patterns: []

Makefile

@@ -156,8 +156,8 @@ bump-version:
 
 jjnew:
 	@echo "$(YELLOW)→ Creating a new commit...$(NC)"
-	@echo "$(BLUE)→ Documenting current commit...$(NC)"
+	@echo "$(BLUE)→ Documenting undocumented commits...$(NC)"
-	@jj auto-describe
+	@jj auto-doc
 	@echo "$(BLUE)→ Done.$(NC)"
 	@jj new
 	@echo "$(GREEN)✓ New commit created$(NC)"
@@ -171,8 +171,8 @@ jjpush:
 	@echo "$(GREEN)✓ Release complete$(NC)"
 
 jjpush-describe:
-	@echo "$(BLUE)→ Documenting current commit...$(NC)"
+	@echo "$(BLUE)→ Documenting undocumented commits...$(NC)"
-	@jj auto-describe
+	@jj auto-doc
 
 jjpush-bump:
 	@echo "$(BLUE)→ Creating new commit for version bump...$(NC)"

cmd/obitools/obipairing/main.go

@@ -37,6 +37,11 @@ func main() {
 
 	optionParser(os.Args)
 
+	if !obipairing.CLIHasPairedFiles() {
+		log.Error("You must provide both a forward file (-F) and a reverse file (-R)")
+		os.Exit(1)
+	}
+
 	obidefault.SetStrictReadWorker(2)
 	obidefault.SetStrictWriteWorker(2)
 	pairs, err := obipairing.CLIPairedSequence()

cmd/obitools/obitagpcr/main.go

@@ -1,6 +1,7 @@
 package main
 
 import (
+	"fmt"
 	"os"
 
 	log "github.com/sirupsen/logrus"
@@ -8,6 +9,7 @@ import (
 	"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/obimultiplex"
 	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obipairing"
 	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obitagpcr"
 	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiutils"
@@ -39,6 +41,17 @@ func main() {
 		obitagpcr.OptionSet)
 
 	optionParser(os.Args)
+
+	if obimultiplex.CLIAskConfigTemplate() {
+		fmt.Print(obimultiplex.CLIConfigTemplate())
+		os.Exit(0)
+	}
+
+	if !obipairing.CLIHasPairedFiles() {
+		log.Error("You must provide both a forward file (-F) and a reverse file (-R)")
+		os.Exit(1)
+	}
+
 	pairs, err := obipairing.CLIPairedSequence()
 
 	if err != nil {

entities.json

@@ -0,0 +1,10 @@
+{
+  "people": [
+    "Software",
+    "Agreement",
+    "Module"
+  ],
+  "projects": [
+    "Code"
+  ]
+}

go.mod

@@ -6,7 +6,7 @@ require (
 	github.com/DavidGamba/go-getoptions v0.33.0
 	github.com/PaesslerAG/gval v1.2.4
 	github.com/barkimedes/go-deepcopy v0.0.0-20220514131651-17c30cfc62df
-	github.com/buger/jsonparser v1.1.1
+	github.com/buger/jsonparser v1.1.2
 	github.com/chen3feng/stl4go v0.1.1
 	github.com/dlclark/regexp2 v1.11.5
 	github.com/goccy/go-json v0.10.6

go.sum

@@ -6,8 +6,8 @@ github.com/PaesslerAG/jsonpath v0.1.0 h1:gADYeifvlqK3R3i2cR5B4DGgxLXIPb3TRTH1mGi
 github.com/PaesslerAG/jsonpath v0.1.0/go.mod h1:4BzmtoM/PI8fPO4aQGIusjGxGir2BzcV0grWtFzq1Y8=
 github.com/barkimedes/go-deepcopy v0.0.0-20220514131651-17c30cfc62df h1:GSoSVRLoBaFpOOds6QyY1L8AX7uoY+Ln3BHc22W40X0=
 github.com/barkimedes/go-deepcopy v0.0.0-20220514131651-17c30cfc62df/go.mod h1:hiVxq5OP2bUGBRNS3Z/bt/reCLFNbdcST6gISi1fiOM=
-github.com/buger/jsonparser v1.1.1 h1:2PnMjfWD7wBILjqQbt530v576A/cAbQvEW9gGIpYMUs=
+github.com/buger/jsonparser v1.1.2 h1:frqHqw7otoVbk5M8LlE/L7HTnIq2v9RX6EJ48i9AxJk=
-github.com/buger/jsonparser v1.1.1/go.mod h1:6RYKKt7H4d4+iWqouImQ9R2FZql3VbhNgx27UK13J/0=
+github.com/buger/jsonparser v1.1.2/go.mod h1:6RYKKt7H4d4+iWqouImQ9R2FZql3VbhNgx27UK13J/0=
 github.com/chen3feng/stl4go v0.1.1 h1:0L1+mDw7pomftKDruM23f1mA7miavOj6C6MZeadzN2Q=
 github.com/chen3feng/stl4go v0.1.1/go.mod h1:5ml3psLgETJjRJnMbPE+JiHLrCpt+Ajc2weeTECXzWU=
 github.com/chengxilo/virtualterm v1.0.4 h1:Z6IpERbRVlfB8WkOmtbHiDbBANU7cimRIof7mk9/PwM=

obitests/obitools/obiconvert/out_ecotag.fasta

@@ -0,0 +1,24 @@
+>HELIUM_000100422_612GNAAXX:7:118:3572:14633#0/1_sub[28..126] {"count":10172,"merged_sample":{"26a_F040644":10172},"obitag_bestid":0.9797979797979798,"obitag_bestmatch":"AY227529","obitag_match_count":1,"obitag_rank":"genus","obitag_similarity_method":"lcs","taxid":"taxon:9992 [Marmota]@genus"}
+ttagccctaaacataaacattcaataaacaagaatgttcgccagagtactactagcaaca
+gcctgaaactcaaaggacttggcggtgctttacatccct
+>HELIUM_000100422_612GNAAXX:7:99:9351:13090#0/1_sub[28..127] {"count":260,"merged_sample":{"29a_F260619":260},"obitag_bestid":0.9405940594059405,"obitag_bestmatch":"AF154263","obitag_match_count":9,"obitag_rank":"infraorder","obitag_similarity_method":"lcs","taxid":"taxon:35500 [Pecora]@infraorder"}
+ttagccctaaacacaaataattacacaaacaaaattgttcaccagagtactagcggcaac
+agcttaaaactcaaaggacttggcggtgctttataccctt
+>HELIUM_000100422_612GNAAXX:7:108:10111:9078#0/1_sub[28..127] {"count":7146,"merged_sample":{"13a_F730603":7146},"obitag_bestid":1,"obitag_bestmatch":"AB245427","obitag_match_count":1,"obitag_rank":"species","obitag_similarity_method":"lcs","taxid":"taxon:9860 [Cervus elaphus]@species"}
+ctagccttaaacacaaatagttatgcaaacaaaactattcgccagagtactaccggcaat
+agcttaaaactcaaaggacttggcggtgctttataccctt
+>HELIUM_000100422_612GNAAXX:7:38:14204:12725#0/1_sub[28..126] {"count":87,"merged_sample":{"26a_F040644":87},"obitag_bestid":0.9494949494949495,"obitag_bestmatch":"AY227530","obitag_match_count":2,"obitag_rank":"tribe","obitag_similarity_method":"lcs","taxid":"taxon:337730 [Marmotini]@tribe"}
+ttagccctaaacataaacattcaataaacaagaatgttcgccagaggactactagcaata
+gcttaaaactcaaaggacttggcggtgctttatatccct
+>HELIUM_000100422_612GNAAXX:7:30:9942:4495#0/1_sub[28..126] {"count":95,"merged_sample":{"26a_F040644":11,"29a_F260619":84},"obitag_bestid":0.9595959595959596,"obitag_bestmatch":"AC187326","obitag_match_count":1,"obitag_rank":"subspecies","obitag_similarity_method":"lcs","taxid":"taxon:9615 [Canis lupus familiaris]@subspecies"}
+ttagccctaaacataagctattccataacaaaataattcgccagagaactactagcaaca
+gattaaacctcaaaggacttggcagtgctttatacccct
+>HELIUM_000100422_612GNAAXX:7:51:16702:19393#0/1_sub[28..127] {"count":12004,"merged_sample":{"15a_F730814":7465,"29a_F260619":4539},"obitag_bestid":1,"obitag_bestmatch":"AJ885202","obitag_match_count":1,"obitag_rank":"species","obitag_similarity_method":"lcs","taxid":"taxon:9858 [Capreolus capreolus]@species"}
+ttagccctaaacacaagtaattaatataacaaaattattcgccagagtactaccggcaat
+agcttaaaactcaaaggacttggcggtgctttataccctt
+>HELIUM_000100422_612GNAAXX:7:84:14502:1617#0/1_sub[28..127] {"count":319,"merged_sample":{"29a_F260619":319},"obitag_bestid":1,"obitag_bestmatch":"AJ972683","obitag_match_count":1,"obitag_rank":"species","obitag_similarity_method":"lcs","taxid":"taxon:9858 [Capreolus capreolus]@species"}
+ttagccctaaacacaagtaattattataacaaaattattcgccagagtactaccggcaat
+agcttaaaactcaaaggacttggcggtgctttataccctt
+>HELIUM_000100422_612GNAAXX:7:50:10637:6527#0/1_sub[28..126] {"count":366,"merged_sample":{"13a_F730603":13,"15a_F730814":5,"26a_F040644":347,"29a_F260619":1},"obitag_bestid":1,"obitag_bestmatch":"AB048590","obitag_match_count":1,"obitag_rank":"genus","obitag_similarity_method":"lcs","taxid":"taxon:9611 [Canis]@genus"}
+ttagccctaaacatagataattttacaacaaaataattcgccagaggactactagcaata
+gcttaaaactcaaaggacttggcggtgctttatatccct

obitests/obitools/obiconvert/taxonomy.csv

@@ -0,0 +1,48 @@
+taxid,parent,taxonomic_rank,scientific_name
+taxon:1 [root]@no rank,taxon:1 [root]@no rank,no rank,root
+taxon:131567 [cellular organisms]@cellular root,taxon:1 [root]@no rank,cellular root,cellular organisms
+taxon:2759 [Eukaryota]@domain,taxon:131567 [cellular organisms]@cellular root,domain,Eukaryota
+taxon:33154 [Opisthokonta]@clade,taxon:2759 [Eukaryota]@domain,clade,Opisthokonta
+taxon:33208 [Metazoa]@kingdom,taxon:33154 [Opisthokonta]@clade,kingdom,Metazoa
+taxon:6072 [Eumetazoa]@clade,taxon:33208 [Metazoa]@kingdom,clade,Eumetazoa
+taxon:33213 [Bilateria]@clade,taxon:6072 [Eumetazoa]@clade,clade,Bilateria
+taxon:33511 [Deuterostomia]@clade,taxon:33213 [Bilateria]@clade,clade,Deuterostomia
+taxon:7711 [Chordata]@phylum,taxon:33511 [Deuterostomia]@clade,phylum,Chordata
+taxon:89593 [Craniata]@subphylum,taxon:7711 [Chordata]@phylum,subphylum,Craniata
+taxon:7742 [Vertebrata]@clade,taxon:89593 [Craniata]@subphylum,clade,Vertebrata
+taxon:7776 [Gnathostomata]@clade,taxon:7742 [Vertebrata]@clade,clade,Gnathostomata
+taxon:117570 [Teleostomi]@clade,taxon:7776 [Gnathostomata]@clade,clade,Teleostomi
+taxon:117571 [Euteleostomi]@clade,taxon:117570 [Teleostomi]@clade,clade,Euteleostomi
+taxon:8287 [Sarcopterygii]@superclass,taxon:117571 [Euteleostomi]@clade,superclass,Sarcopterygii
+taxon:1338369 [Dipnotetrapodomorpha]@clade,taxon:8287 [Sarcopterygii]@superclass,clade,Dipnotetrapodomorpha
+taxon:32523 [Tetrapoda]@clade,taxon:1338369 [Dipnotetrapodomorpha]@clade,clade,Tetrapoda
+taxon:32524 [Amniota]@clade,taxon:32523 [Tetrapoda]@clade,clade,Amniota
+taxon:40674 [Mammalia]@class,taxon:32524 [Amniota]@clade,class,Mammalia
+taxon:32525 [Theria]@clade,taxon:40674 [Mammalia]@class,clade,Theria
+taxon:9347 [Eutheria]@clade,taxon:32525 [Theria]@clade,clade,Eutheria
+taxon:1437010 [Boreoeutheria]@clade,taxon:9347 [Eutheria]@clade,clade,Boreoeutheria
+taxon:314146 [Euarchontoglires]@superorder,taxon:1437010 [Boreoeutheria]@clade,superorder,Euarchontoglires
+taxon:314145 [Laurasiatheria]@superorder,taxon:1437010 [Boreoeutheria]@clade,superorder,Laurasiatheria
+taxon:33554 [Carnivora]@order,taxon:314145 [Laurasiatheria]@superorder,order,Carnivora
+taxon:91561 [Artiodactyla]@order,taxon:314145 [Laurasiatheria]@superorder,order,Artiodactyla
+taxon:314147 [Glires]@clade,taxon:314146 [Euarchontoglires]@superorder,clade,Glires
+taxon:9845 [Ruminantia]@suborder,taxon:91561 [Artiodactyla]@order,suborder,Ruminantia
+taxon:35500 [Pecora]@infraorder,taxon:9845 [Ruminantia]@suborder,infraorder,Pecora
+taxon:9989 [Rodentia]@order,taxon:314147 [Glires]@clade,order,Rodentia
+taxon:379584 [Caniformia]@suborder,taxon:33554 [Carnivora]@order,suborder,Caniformia
+taxon:9608 [Canidae]@family,taxon:379584 [Caniformia]@suborder,family,Canidae
+taxon:9850 [Cervidae]@family,taxon:35500 [Pecora]@infraorder,family,Cervidae
+taxon:9881 [Odocoileinae]@subfamily,taxon:9850 [Cervidae]@family,subfamily,Odocoileinae
+taxon:33553 [Sciuromorpha]@suborder,taxon:9989 [Rodentia]@order,suborder,Sciuromorpha
+taxon:55153 [Sciuridae]@family,taxon:33553 [Sciuromorpha]@suborder,family,Sciuridae
+taxon:34878 [Cervinae]@subfamily,taxon:9850 [Cervidae]@family,subfamily,Cervinae
+taxon:9611 [Canis]@genus,taxon:9608 [Canidae]@family,genus,Canis
+taxon:9857 [Capreolus]@genus,taxon:9881 [Odocoileinae]@subfamily,genus,Capreolus
+taxon:9612 [Canis lupus]@species,taxon:9611 [Canis]@genus,species,Canis lupus
+taxon:337726 [Xerinae]@subfamily,taxon:55153 [Sciuridae]@family,subfamily,Xerinae
+taxon:9859 [Cervus]@genus,taxon:34878 [Cervinae]@subfamily,genus,Cervus
+taxon:337730 [Marmotini]@tribe,taxon:337726 [Xerinae]@subfamily,tribe,Marmotini
+taxon:9992 [Marmota]@genus,taxon:337730 [Marmotini]@tribe,genus,Marmota
+taxon:9860 [Cervus elaphus]@species,taxon:9859 [Cervus]@genus,species,Cervus elaphus
+taxon:9615 [Canis lupus familiaris]@subspecies,taxon:9612 [Canis lupus]@species,subspecies,Canis lupus familiaris
+taxon:9858 [Capreolus capreolus]@species,taxon:9857 [Capreolus]@genus,species,Capreolus capreolus

obitests/obitools/obiconvert/test.sh

@@ -134,6 +134,130 @@ else
     ((failed++))
 fi
 
+
+# ------------------------------------------------------------------
+# --raw-taxid tests (no taxonomy loaded)
+# ------------------------------------------------------------------
+
+# Running test
+((ntest++))
+if obiconvert --raw-taxid "${TEST_DIR}/out_ecotag.fasta" \
+              > "${TMPDIR}/raw_taxid.fasta" 2>/dev/null
+then
+    log "$MCMD --raw-taxid: running OK"
+    ((success++))
+else
+    log "$MCMD --raw-taxid: running failed"
+    ((failed++))
+fi
+
+# Taxids must be bare numbers — no full-format "taxon:ID [Name]@rank" strings
+((ntest++))
+if grep '"taxid"' "${TMPDIR}/raw_taxid.fasta" | grep -qv '"taxid":"[0-9][0-9]*"'
+then
+    log "$MCMD --raw-taxid: taxid format check failed (full-format taxid found)"
+    ((failed++))
+else
+    log "$MCMD --raw-taxid: taxid format OK (all taxids are bare numbers)"
+    ((success++))
+fi
+
+# --raw-taxid is idempotent: piping through a second obiconvert --raw-taxid must
+# produce bit-for-bit identical output.
+((ntest++))
+if obiconvert --raw-taxid "${TMPDIR}/raw_taxid.fasta" \
+              > "${TMPDIR}/raw_taxid2.fasta" 2>/dev/null
+then
+    log "$MCMD --raw-taxid piped: running OK"
+    ((success++))
+else
+    log "$MCMD --raw-taxid piped: running failed"
+    ((failed++))
+fi
+
+((ntest++))
+if diff "${TMPDIR}/raw_taxid.fasta" \
+        "${TMPDIR}/raw_taxid2.fasta" > /dev/null
+then
+    log "$MCMD --raw-taxid piped: idempotency OK"
+    ((success++))
+else
+    log "$MCMD --raw-taxid piped: idempotency failed (outputs differ)"
+    ((failed++))
+fi
+
+
+# ------------------------------------------------------------------
+# --taxonomy tests (full-format taxid, no --raw-taxid)
+# ------------------------------------------------------------------
+
+# Running test
+((ntest++))
+if obiconvert --taxonomy "${TEST_DIR}/taxonomy.csv" \
+              "${TEST_DIR}/out_ecotag.fasta" \
+              > "${TMPDIR}/taxo.fasta" 2>/dev/null
+then
+    log "$MCMD --taxonomy: running OK"
+    ((success++))
+else
+    log "$MCMD --taxonomy: running failed"
+    ((failed++))
+fi
+
+# Taxids must be in full "taxon:ID [Name]@rank" format
+((ntest++))
+if grep '"taxid"' "${TMPDIR}/taxo.fasta" | grep -q '"taxid":"taxon:[0-9]'
+then
+    log "$MCMD --taxonomy: taxid format OK (full-format taxids present)"
+    ((success++))
+else
+    log "$MCMD --taxonomy: taxid format check failed (no full-format taxid found)"
+    ((failed++))
+fi
+
+
+# ------------------------------------------------------------------
+# --raw-taxid --taxonomy tests
+# ------------------------------------------------------------------
+
+# Running test
+((ntest++))
+if obiconvert --raw-taxid --taxonomy "${TEST_DIR}/taxonomy.csv" \
+              "${TEST_DIR}/out_ecotag.fasta" \
+              > "${TMPDIR}/raw_taxid_taxo.fasta" 2>/dev/null
+then
+    log "$MCMD --raw-taxid --taxonomy: running OK"
+    ((success++))
+else
+    log "$MCMD --raw-taxid --taxonomy: running failed"
+    ((failed++))
+fi
+
+# Taxids must be bare numbers even when taxonomy is loaded
+((ntest++))
+if grep '"taxid"' "${TMPDIR}/raw_taxid_taxo.fasta" | grep -qv '"taxid":"[0-9][0-9]*"'
+then
+    log "$MCMD --raw-taxid --taxonomy: taxid format check failed (full-format taxid found)"
+    ((failed++))
+else
+    log "$MCMD --raw-taxid --taxonomy: taxid format OK (all taxids are bare numbers)"
+    ((success++))
+fi
+
+# --raw-taxid with or without taxonomy must yield identical taxid values
+((ntest++))
+if diff <(grep '"taxid"' "${TMPDIR}/raw_taxid.fasta" | grep -o '"taxid":"[^"]*"' | sort) \
+        <(grep '"taxid"' "${TMPDIR}/raw_taxid_taxo.fasta" | grep -o '"taxid":"[^"]*"' | sort) \
+        > /dev/null
+then
+    log "$MCMD --raw-taxid vs --raw-taxid --taxonomy: taxid values match OK"
+    ((success++))
+else
+    log "$MCMD --raw-taxid vs --raw-taxid --taxonomy: taxid values differ (unexpected)"
+    ((failed++))
+fi
+
+
 #########################################
 #
 # At the end of the tests

obitests/obitools/obiconvert/xxx

@@ -0,0 +1,24 @@
+>HELIUM_000100422_612GNAAXX:7:118:3572:14633#0/1_sub[28..126] {"count":10172,"merged_sample":{"26a_F040644":10172},"obitag_bestid":0.9797979797979798,"obitag_bestmatch":"AY227529","obitag_match_count":1,"obitag_rank":"genus","obitag_similarity_method":"lcs","taxid":"taxon:9992 [Marmota]@genus"}
+ttagccctaaacataaacattcaataaacaagaatgttcgccagagtactactagcaaca
+gcctgaaactcaaaggacttggcggtgctttacatccct
+>HELIUM_000100422_612GNAAXX:7:99:9351:13090#0/1_sub[28..127] {"count":260,"merged_sample":{"29a_F260619":260},"obitag_bestid":0.9405940594059405,"obitag_bestmatch":"AF154263","obitag_match_count":9,"obitag_rank":"infraorder","obitag_similarity_method":"lcs","taxid":"taxon:35500 [Pecora]@infraorder"}
+ttagccctaaacacaaataattacacaaacaaaattgttcaccagagtactagcggcaac
+agcttaaaactcaaaggacttggcggtgctttataccctt
+>HELIUM_000100422_612GNAAXX:7:108:10111:9078#0/1_sub[28..127] {"count":7146,"merged_sample":{"13a_F730603":7146},"obitag_bestid":1,"obitag_bestmatch":"AB245427","obitag_match_count":1,"obitag_rank":"species","obitag_similarity_method":"lcs","taxid":"taxon:9860 [Cervus elaphus]@species"}
+ctagccttaaacacaaatagttatgcaaacaaaactattcgccagagtactaccggcaat
+agcttaaaactcaaaggacttggcggtgctttataccctt
+>HELIUM_000100422_612GNAAXX:7:38:14204:12725#0/1_sub[28..126] {"count":87,"merged_sample":{"26a_F040644":87},"obitag_bestid":0.9494949494949495,"obitag_bestmatch":"AY227530","obitag_match_count":2,"obitag_rank":"tribe","obitag_similarity_method":"lcs","taxid":"taxon:337730 [Marmotini]@tribe"}
+ttagccctaaacataaacattcaataaacaagaatgttcgccagaggactactagcaata
+gcttaaaactcaaaggacttggcggtgctttatatccct
+>HELIUM_000100422_612GNAAXX:7:30:9942:4495#0/1_sub[28..126] {"count":95,"merged_sample":{"26a_F040644":11,"29a_F260619":84},"obitag_bestid":0.9595959595959596,"obitag_bestmatch":"AC187326","obitag_match_count":1,"obitag_rank":"subspecies","obitag_similarity_method":"lcs","taxid":"taxon:9615 [Canis lupus familiaris]@subspecies"}
+ttagccctaaacataagctattccataacaaaataattcgccagagaactactagcaaca
+gattaaacctcaaaggacttggcagtgctttatacccct
+>HELIUM_000100422_612GNAAXX:7:51:16702:19393#0/1_sub[28..127] {"count":12004,"merged_sample":{"15a_F730814":7465,"29a_F260619":4539},"obitag_bestid":1,"obitag_bestmatch":"AJ885202","obitag_match_count":1,"obitag_rank":"species","obitag_similarity_method":"lcs","taxid":"taxon:9858 [Capreolus capreolus]@species"}
+ttagccctaaacacaagtaattaatataacaaaattattcgccagagtactaccggcaat
+agcttaaaactcaaaggacttggcggtgctttataccctt
+>HELIUM_000100422_612GNAAXX:7:84:14502:1617#0/1_sub[28..127] {"count":319,"merged_sample":{"29a_F260619":319},"obitag_bestid":1,"obitag_bestmatch":"AJ972683","obitag_match_count":1,"obitag_rank":"species","obitag_similarity_method":"lcs","taxid":"taxon:9858 [Capreolus capreolus]@species"}
+ttagccctaaacacaagtaattattataacaaaattattcgccagagtactaccggcaat
+agcttaaaactcaaaggacttggcggtgctttataccctt
+>HELIUM_000100422_612GNAAXX:7:50:10637:6527#0/1_sub[28..126] {"count":366,"merged_sample":{"13a_F730603":13,"15a_F730814":5,"26a_F040644":347,"29a_F260619":1},"obitag_bestid":1,"obitag_bestmatch":"AB048590","obitag_match_count":1,"obitag_rank":"genus","obitag_similarity_method":"lcs","taxid":"taxon:9611 [Canis]@genus"}
+ttagccctaaacatagataattttacaacaaaataattcgccagaggactactagcaata
+gcttaaaactcaaaggacttggcggtgctttatatccct

pkg/obialign/locatepattern.go

@@ -28,10 +28,18 @@ func buffIndex(i, j, width int) int {
 //
 // The function returns the start and end positions of the best
 // match, as well as the number of errors in the best match.
+//
+// When the sequence is too short relative to the pattern for the
+// backtracking to reconstruct a valid alignment (e.g. the pattern
+// is longer than the sequence, or the match sits too close to a
+// sequence boundary), no reliable position can be computed. In that
+// case the function returns the sentinel (-1, -1, -1) instead of a
+// guessed, potentially wrong, position: callers must treat this as
+// "no match" rather than use the returned coordinates.
 func LocatePattern(id string, pattern, sequence []byte) (int, int, int) {
 
-	if len(pattern) >= len(sequence) {
+	if len(sequence) == 0 {
-		log.Panicf("Sequence %s:Pattern %s must be shorter than sequence %s", id, pattern, sequence)
+		log.Panicf("Sequence %s:Pattern %s must not be empty", id, pattern)
 	}
 
 	// Pattern spreads over the columns
@@ -158,5 +166,18 @@ func LocatePattern(id string, pattern, sequence []byte) (int, int, int) {
 	// obilog.Warnf("from : %d to: %d error: %d match: %v",
 	// 	i, end+1, -buffer[buffIndex(len(sequence)-1, len(pattern)-1, width)],
 	// 	string(sequence[i:(end+1)]))
+
+	if i < 0 || end == -1 {
+		// i < 0: the backtracking ran off the start of the sequence
+		// without fully consuming the pattern.
+		// end == -1: the backtracking loop never ran at all (e.g. a
+		// single-base pattern, jmax == 0), so no alignment boundary
+		// was ever established.
+		// Either way, no valid alignment exists for this (pattern,
+		// sequence) pair: signal it explicitly instead of returning
+		// an out-of-bounds or uncomputed position.
+		return -1, -1, -1
+	}
+
 	return i, end + 1, -buffer[buffIndex(len(sequence)-1, len(pattern)-1, width)]
 }

pkg/obialign/locatepattern_test.go

@@ -0,0 +1,123 @@
+package obialign
+
+import (
+	"math/rand"
+	"testing"
+)
+
+func TestLocatePatternNormal(t *testing.T) {
+	// Pattern fully and exactly present in the middle of a longer sequence.
+	start, end, nerr := LocatePattern("id", []byte("ACGT"), []byte("TTTTACGTTTTT"))
+	if start != 4 || end != 8 || nerr != 0 {
+		t.Errorf("got start=%d end=%d nerr=%d, want start=4 end=8 nerr=0", start, end, nerr)
+	}
+}
+
+func TestLocatePatternOneMismatch(t *testing.T) {
+	start, end, nerr := LocatePattern("id", []byte("ACGT"), []byte("TTTTACTTTTTT"))
+	if nerr != 1 {
+		t.Errorf("got nerr=%d, want 1 (start=%d end=%d)", nerr, start, end)
+	}
+}
+
+// The real-world case that used to panic: pattern longer than the sequence
+// fragment extracted for indel relocation.
+func TestLocatePatternPatternLongerThanSequence(t *testing.T) {
+	start, end, nerr := LocatePattern("id",
+		[]byte("GGGCAATCCTGAGCCAAATC"),
+		[]byte("tcctgagccaaatcacgtt"))
+
+	if start != -1 || end != -1 || nerr != -1 {
+		t.Errorf("got start=%d end=%d nerr=%d, want the (-1,-1,-1) sentinel", start, end, nerr)
+	}
+}
+
+func TestLocatePatternSequenceLengthOne(t *testing.T) {
+	start, end, nerr := LocatePattern("id", []byte("AB"), []byte("A"))
+
+	if start < 0 || end < 0 || nerr < 0 {
+		t.Fatalf("got start=%d end=%d nerr=%d, expected a valid (non-sentinel) result", start, end, nerr)
+	}
+	if start != 0 || end != 1 || nerr != 1 {
+		t.Errorf("got start=%d end=%d nerr=%d, want start=0 end=1 nerr=1", start, end, nerr)
+	}
+}
+
+// A pattern much longer than the sequence can still yield a mathematically
+// valid (in-bounds) alignment: the extra pattern length is absorbed as gaps,
+// driving the error count high enough that the caller's maxerr threshold
+// rejects it. The function itself must still return consistent bounds.
+func TestLocatePatternPatternMuchLongerThanSequence(t *testing.T) {
+	start, end, nerr := LocatePattern("id", []byte("ACGTACGTACGTACGTACGT"), []byte("ACG"))
+
+	isSentinel := start == -1 && end == -1 && nerr == -1
+	isValid := start >= 0 && end > start && end <= 3 && nerr >= 0
+
+	if !isSentinel && !isValid {
+		t.Errorf("got start=%d end=%d nerr=%d, want either the sentinel or consistent in-bounds values", start, end, nerr)
+	}
+}
+
+func TestLocatePatternNeverReturnsOutOfBounds(t *testing.T) {
+	patterns := []string{"A", "AC", "ACG", "ACGT", "ACGTA", "ACGTAC", "ACGTACG", "ACGTACGT"}
+	sequences := []string{"A", "AC", "ACG", "ACGT", "ACGTA", "ACGTAC", "ACGTACG", "ACGTACGT"}
+
+	for _, p := range patterns {
+		for _, s := range sequences {
+			start, end, nerr := LocatePattern("id", []byte(p), []byte(s))
+
+			if start == -1 && end == -1 && nerr == -1 {
+				// Explicit "no reliable match" sentinel: always acceptable.
+				continue
+			}
+
+			if start < 0 || end < 0 || start >= end || end > len(s) || nerr < 0 {
+				t.Errorf("pattern=%q sequence=%q -> start=%d end=%d nerr=%d is out of bounds / inconsistent",
+					p, s, start, end, nerr)
+			}
+		}
+	}
+}
+
+// Randomized property test over a wide range of pattern/sequence length
+// combinations, including pattern >= sequence, to make sure the function
+// never panics and never returns anything but the sentinel or fully
+// consistent, in-bounds coordinates.
+func TestLocatePatternRandomizedNeverInvalid(t *testing.T) {
+	const bases = "ACGT"
+	rng := rand.New(rand.NewSource(42))
+
+	randSeq := func(n int) []byte {
+		b := make([]byte, n)
+		for i := range b {
+			b[i] = bases[rng.Intn(len(bases))]
+		}
+		return b
+	}
+
+	for trial := 0; trial < 5000; trial++ {
+		patLen := 1 + rng.Intn(15)
+		seqLen := 1 + rng.Intn(15)
+
+		pattern := randSeq(patLen)
+		sequence := randSeq(seqLen)
+
+		func() {
+			defer func() {
+				if r := recover(); r != nil {
+					t.Fatalf("panic for pattern=%q sequence=%q: %v", pattern, sequence, r)
+				}
+			}()
+
+			start, end, nerr := LocatePattern("id", pattern, sequence)
+
+			isSentinel := start == -1 && end == -1 && nerr == -1
+			isValid := start >= 0 && end > start && end <= len(sequence) && nerr >= 0
+
+			if !isSentinel && !isValid {
+				t.Errorf("pattern=%q sequence=%q -> start=%d end=%d nerr=%d is neither the sentinel nor consistent",
+					pattern, sequence, start, end, nerr)
+			}
+		}()
+	}
+}

pkg/obiapat/pattern.go

@@ -373,6 +373,7 @@ func (pattern ApatPattern) BestMatch(sequence ApatSequence, begin, length int) (
 
 	cpattern := (*[1 << 30]byte)(unsafe.Pointer(pattern.pointer.pointer.cpat))
 	frg := sequence.pointer.reference.Sequence()[start:end]
+	fragStart := start
 
 	log.Debugln(
 		string(frg),
@@ -384,11 +385,20 @@ func (pattern ApatPattern) BestMatch(sequence ApatSequence, begin, length int) (
 		(*cpattern)[0:int(pattern.pointer.pointer.patlen)],
 		frg)
 
-	// olderr := m[2]
+	if from < 0 {
+		// obialign.LocatePattern could not reconstruct a reliable
+		// alignment (e.g. the fragment is too short relative to the
+		// pattern). Reporting a guessed position would risk placing
+		// the primer boundary incorrectly, so treat it as no match
+		// at all rather than falling back to an unrefined position.
+		matched = false
+		log.Debugln("No reliable indel relocation, discarding match", sequence.pointer.reference.Id())
+		return
+	}
 
 	nerr = score
-	start = start + from
+	start = fragStart + from
-	end = start + to
+	end = fragStart + to
 	log.Debugf("BestMatch on %s : score=%d [%d..%d]", sequence.pointer.reference.Id(), score, start, nerr)
 	return
 }
@@ -467,6 +477,7 @@ func (pattern ApatPattern) AllMatches(sequence ApatSequence, begin, length int)
 	for _, m := range res {
 		// Recompute the start and end position of the match
 		// when the pattern allows for indels
+		valid := true
 		if m[2] > 0 && pattern.pointer.pointer.hasIndel {
 			// obilog.Warnf("Locating indel on sequence %s[%s]", sequence.pointer.reference.Id(), pattern.String())
 			start := m[0] - m[2]*2
@@ -485,6 +496,15 @@ func (pattern ApatPattern) AllMatches(sequence ApatSequence, begin, length int)
 				(*cpattern)[0:int(pattern.pointer.pointer.patlen)],
 				frg)
 
+			if pb < 0 {
+				// obialign.LocatePattern could not reconstruct a
+				// reliable alignment (e.g. the match sits too close
+				// to a sequence end for the fragment to be usable).
+				// Reporting a guessed position risks placing the
+				// primer boundary incorrectly, so drop the match
+				// entirely instead of keeping an unrefined guess.
+				valid = false
+			} else {
 				// olderr := m[2]
 				m[2] = score
 				m[0] = start + pb
@@ -493,8 +513,9 @@ func (pattern ApatPattern) AllMatches(sequence ApatSequence, begin, length int)
 				// obilog.Warnf("seq[%d@%d:%d] %d: %s %d - %s:%s:%s", i, m[0], m[1], olderr, sequence.pointer.reference.Id(), score,
 				// 	frg, (*cpattern)[0:int(pattern.pointer.pointer.patlen)], sequence.pointer.reference.Sequence()[m[0]:m[1]])
 			}
+		}
 
-		if int(pattern.pointer.pointer.maxerr) >= m[2] {
+		if valid && int(pattern.pointer.pointer.maxerr) >= m[2] {
 			res[j] = m
 			j++
 		}

pkg/obiformats/fastqseq_read.go

@@ -131,7 +131,7 @@ func _storeSequenceQuality(bytes *bytes.Buffer, out *obiseq.BioSequence, quality
 	out.SetQualities(q)
 }
 
-func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(string, io.Reader) (obiseq.BioSequenceSlice, error) {
+func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool, fileName string) func(string, io.Reader) (obiseq.BioSequenceSlice, error) {
 	parser := func(source string, input io.Reader) (obiseq.BioSequenceSlice, error) {
 
 		var identifier string
@@ -160,12 +160,12 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str
 					// Beginning of sequence
 					state = 1
 				} else {
-					log.Fatalf("%s : sequence entry is not starting with @", source)
+					log.Fatalf("file %s: sequence entry is not starting with @", fileName)
 				}
 			case 1: // Beginning of identifier (Mandatory)
 				if is_sep {
 					// No identifier -> ERROR
-					log.Fatalf("%s : sequence identifier is empty", source)
+					log.Fatalf("file %s: sequence identifier is empty", fileName)
 				} else {
 					// Beginning of identifier
 					state = 2
@@ -221,7 +221,7 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str
 					// End of sequence
 					rawseq := seqBytes.Bytes()
 					if len(rawseq) == 0 {
-						log.Fatalf("@%s[%s] : sequence is empty", identifier, source)
+						log.Fatalf("file %s: record @%s has an empty sequence line", fileName, identifier)
 					}
 					s := obiseq.NewBioSequence(identifier, rawseq, definition)
 					s.SetSource(source)
@@ -241,8 +241,8 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str
 						context = append(
 							append([]byte{previous}, C),
 							context...)
-						log.Fatalf("%s [%s]: sequence contains invalid character %c (%s)",
+						log.Fatalf("file %s: record @%s contains invalid character %c (%s)",
-							source, identifier, C, string(context))
+							fileName, identifier, C, string(context))
 					}
 				}
 			case 7:
@@ -251,7 +251,7 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str
 				} else if C == '+' {
 					state = 8
 				} else {
-					log.Fatalf("@%s[%s] : sequence data not followed by a line starting with + but a %c", identifier, source, C)
+					log.Fatalf("file %s: record @%s: sequence data not followed by a line starting with + but a %c", fileName, identifier, C)
 				}
 			case 8:
 				// State consuming the + internal header line
@@ -282,7 +282,7 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str
 				} else if C == '@' {
 					state = 1
 				} else {
-					log.Fatalf("%s[%s] : sequence record not followed by a line starting with @", identifier, source)
+					log.Fatalf("file %s: record @%s not followed by a line starting with @", fileName, identifier)
 				}
 
 			}
@@ -304,7 +304,7 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str
 }
 
 // FastqChunkParserRope parses a FASTQ chunk directly from a rope without Pack().
-func FastqChunkParserRope(source string, rope *PieceOfChunk, quality_shift byte, with_quality, UtoT bool) (obiseq.BioSequenceSlice, error) {
+func FastqChunkParserRope(source string, rope *PieceOfChunk, quality_shift byte, with_quality, UtoT bool, fileName string) (obiseq.BioSequenceSlice, error) {
 	scanner := newRopeScanner(rope)
 	sequences := obiseq.MakeBioSequenceSlice(100)[:0]
 
@@ -334,7 +334,7 @@ func FastqChunkParserRope(source string, rope *PieceOfChunk, quality_shift byte,
 		// Line 2: sequence
 		sline := scanner.ReadLine()
 		if sline == nil {
-			log.Fatalf("@%s[%s]: unexpected EOF after header", id, source)
+			log.Fatalf("file %s: record @%s is truncated (header line with no sequence line following) — the FASTQ file appears incomplete", fileName, id)
 		}
 		seqDest := make([]byte, len(sline))
 		w := 0
@@ -350,7 +350,7 @@ func FastqChunkParserRope(source string, rope *PieceOfChunk, quality_shift byte,
 		}
 		seqDest = seqDest[:w]
 		if len(seqDest) == 0 {
-			log.Fatalf("@%s[%s]: sequence is empty", id, source)
+			log.Fatalf("file %s: record @%s has an empty sequence line", fileName, id)
 		}
 
 		// Line 3: + (skip)
@@ -382,16 +382,17 @@ func _ParseFastqFile(
 	out obiiter.IBioSequence,
 	quality_shift byte,
 	with_quality, UtoT bool,
+	fileName string,
 ) {
 
-	parser := FastqChunkParser(quality_shift, with_quality, UtoT)
+	parser := FastqChunkParser(quality_shift, with_quality, UtoT, fileName)
 
 	for chunks := range input {
 		var sequences obiseq.BioSequenceSlice
 		var err error
 
 		if chunks.Rope != nil {
-			sequences, err = FastqChunkParserRope(chunks.Source, chunks.Rope, quality_shift, with_quality, UtoT)
+			sequences, err = FastqChunkParserRope(chunks.Source, chunks.Rope, quality_shift, with_quality, UtoT, fileName)
 		} else {
 			sequences, err = parser(chunks.Source, chunks.Raw)
 		}
@@ -423,6 +424,8 @@ func ReadFastq(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, e
 		false,
 	)
 
+	fileName := opt.FileName()
+
 	for i := 0; i < nworker; i++ {
 		out.Add(1)
 		go _ParseFastqFile(
@@ -431,6 +434,7 @@ func ReadFastq(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, e
 			obidefault.ReadQualitiesShift(),
 			opt.ReadQualities(),
 			opt.UtoT(),
+			fileName,
 		)
 	}
 
@@ -456,7 +460,9 @@ func ReadFastq(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, e
 }
 
 func ReadFastqFromFile(filename string, options ...WithOption) (obiiter.IBioSequence, error) {
-	options = append(options, OptionsSource(obiutils.RemoveAllExt((path.Base(filename)))))
+	options = append(options,
+		OptionsSource(obiutils.RemoveAllExt((path.Base(filename)))),
+		OptionsFileName(filename))
 
 	file, err := obiutils.Ropen(filename)
 

pkg/obiformats/fastseq_json_header.go

@@ -199,47 +199,13 @@ func _parse_json_array_interface(str []byte) ([]interface{}, error) {
 	return values, nil
 }
 
-func _parse_json_header_(header string, sequence *obiseq.BioSequence) string {
+// _parse_json_annotation_field parses a single key/value pair coming from a
+// JSON object (either a FASTA/FASTQ inline JSON header, or the "annotations"
+// field of a JSON sequence record) and applies it to the sequence, special
+// casing the well-known OBITools attributes (id, definition, count, taxid,
+// obiclean_*, merged_*).
+func _parse_json_annotation_field(key []byte, value []byte, dataType jsonparser.ValueType, sequence *obiseq.BioSequence) error {
 	annotations := sequence.Annotations()
-	start := -1
-	stop := -1
-	level := 0
-	lh := len(header)
-	inquote := false
-
-	for i := 0; (i < lh) && (stop < 0); i++ {
-		// fmt.Printf("[%d,%d-%d] : %d (%c) (%d,%c)\n", i, start, stop, header[i], header[i], '{', '{')
-		if level == 0 && header[i] == '{' && !inquote {
-			start = i
-		}
-
-		// TODO: escaped double quotes are not considered
-		if start > -1 && header[i] == '"' {
-			inquote = !inquote
-		}
-
-		if header[i] == '{' && !inquote {
-			level++
-		}
-
-		if header[i] == '}' && !inquote {
-			level--
-		}
-
-		if start >= 0 && level == 0 {
-			stop = i
-		}
-
-	}
-
-	if start < 0 || stop < 0 {
-		return header
-	}
-
-	stop++
-
-	jsonparser.ObjectEach(obiutils.UnsafeBytes(header[start:stop]),
-		func(key []byte, value []byte, dataType jsonparser.ValueType, offset int) error {
 	var err error
 
 	skey := obiutils.UnsafeString(key)
@@ -335,6 +301,49 @@ func _parse_json_header_(header string, sequence *obiseq.BioSequence) string {
 	}
 
 	return err
+}
+
+func _parse_json_header_(header string, sequence *obiseq.BioSequence) string {
+	start := -1
+	stop := -1
+	level := 0
+	lh := len(header)
+	inquote := false
+
+	for i := 0; (i < lh) && (stop < 0); i++ {
+		// fmt.Printf("[%d,%d-%d] : %d (%c) (%d,%c)\n", i, start, stop, header[i], header[i], '{', '{')
+		if level == 0 && header[i] == '{' && !inquote {
+			start = i
+		}
+
+		// TODO: escaped double quotes are not considered
+		if start > -1 && header[i] == '"' {
+			inquote = !inquote
+		}
+
+		if header[i] == '{' && !inquote {
+			level++
+		}
+
+		if header[i] == '}' && !inquote {
+			level--
+		}
+
+		if start >= 0 && level == 0 {
+			stop = i
+		}
+
+	}
+
+	if start < 0 || stop < 0 {
+		return header
+	}
+
+	stop++
+
+	jsonparser.ObjectEach(obiutils.UnsafeBytes(header[start:stop]),
+		func(key []byte, value []byte, dataType jsonparser.ValueType, offset int) error {
+			return _parse_json_annotation_field(key, value, dataType, sequence)
 		},
 	)
 

pkg/obiformats/fastseq_obi_header.go

@@ -17,7 +17,7 @@ import (
 )
 
 var __obi_header_value_string_pattern__ = regexp.MustCompile(`^'\s*([^']*'|"[^"]*")\s*;`)
-var __obi_header_value_numeric_pattern__ = regexp.MustCompile(`^\s*([+-]?\.\d+|[+-]?\d+(\.\d*)?([eE][+-]?\d+)?)\s*;`)
+var __obi_header_value_numeric_pattern__ = regexp.MustCompile(`^\s*[+-]?(\d+(\.\d*)?|\.\d+)([eE][+-]?\d+)?\s*;`)
 var __obi_header_map_int_key__ = regexp.MustCompile("([{,])([0-9]+):")
 
 func __match__dict__(text []byte) []int {
@@ -146,6 +146,65 @@ func __match__key__(text []byte) []int {
 	return []int{} // Not a key
 }
 
+func __match__array__(text []byte) []int {
+
+	state := 0
+	level := 0
+	start := 0
+	instring := byte(0)
+
+	for i, r := range text {
+		if state == 2 {
+			if r == ';' {
+				return []int{start, i + 1}
+			}
+			if r != ' ' && r != '\t' {
+				return []int{}
+			}
+		}
+
+		if state == 0 {
+			if r == '[' {
+				level++
+				state++
+				start = i
+				continue
+			}
+			if r != ' ' && r != '\t' {
+				return []int{}
+			}
+			continue
+		}
+
+		// state == 1: inside the array
+		if instring != 0 {
+			if r == instring {
+				instring = 0
+			}
+			continue
+		}
+
+		if r == '"' || r == '\'' {
+			instring = r
+			continue
+		}
+
+		if r == '[' || r == '{' {
+			level++
+			continue
+		}
+
+		if r == ']' || r == '}' {
+			level--
+			if level == 0 {
+				state++
+			}
+		}
+	}
+
+	return []int{}
+}
+
 func __match__general__(text []byte) []int {
 
 	for i, r := range text {
@@ -242,6 +301,21 @@ func ParseOBIFeatures(text string, annotations obiseq.Annotation) string {
 					stop = m[1] + 1
 				} else {
 
+					// array value
+					m = __match__array__(part)
+					if len(m) > 0 {
+						bvalue = bytes.TrimSpace(part[m[0]:(m[1] - 1)])
+						j := bytes.ReplaceAll(bvalue, []byte("'"), []byte(`"`))
+						j = __obi_header_map_int_key__.ReplaceAll(j, []byte(`$1"$2":`))
+						arr, err := _parse_json_array_interface(j)
+						if err != nil {
+							value = string(bvalue)
+						} else {
+							value = arr
+						}
+						stop = m[1] + 1
+					} else {
+
 						// Generic value
 
 						// m = __obi_header_value_general_pattern__.FindIndex(part)
@@ -264,6 +338,7 @@ func ParseOBIFeatures(text string, annotations obiseq.Annotation) string {
 							// no value
 							break
 						} // End of No value
+					} // End of not array
 				} // End of not dict
 			} // End of not string
 		} // End of not numeric
@@ -272,6 +347,8 @@ func ParseOBIFeatures(text string, annotations obiseq.Annotation) string {
 		case float64:
 			if vt == math.Floor(vt) {
 				annotations[key] = int(vt)
+			} else {
+				annotations[key] = vt
 			}
 		default:
 			annotations[key] = value
@@ -327,9 +404,8 @@ func WriteFastSeqOBIHeade(buffer *bytes.Buffer, sequence *obiseq.BioSequence) {
 					buffer.WriteString(fmt.Sprintf("%s=", key))
 					buffer.Write(tv)
 					buffer.WriteString("; ")
-				case map[string]int,
+				default:
-					map[string]string,
+					if obiutils.IsAMap(value) || obiutils.IsASlice(value) || obiutils.IsAnArray(value) {
-					map[string]interface{}:
 						tv, err := obiutils.JsonMarshal(t)
 						if err != nil {
 							log.Fatalf("Cannot convert %v value", value)
@@ -338,11 +414,12 @@ func WriteFastSeqOBIHeade(buffer *bytes.Buffer, sequence *obiseq.BioSequence) {
 						buffer.WriteString(fmt.Sprintf("%s=", key))
 						buffer.Write(tv)
 						buffer.WriteString("; ")
-				default:
+					} else {
 						buffer.WriteString(fmt.Sprintf("%s=%v; ", key, value))
 					}
 				}
 			}
+		}
 
 		if sequence.HasDefinition() {
 			buffer.WriteByte(' ')

pkg/obiformats/fastseq_write_fasta.go

@@ -90,6 +90,9 @@ func FormatFastaBatch(batch obiiter.BioSequenceBatch, formater FormatHeader, ski
 	log.Debugf("FormatFastaBatch: #%d : %d seqs", batch.Order(), batch.Len())
 
 	for _, seq := range batch.Slice() {
+		if len(seq.Id()) == 0 {
+			log.Fatalf("Sequence identifier is empty")
+		}
 		if seq.Len() > 0 {
 			// Write header directly into bs — no intermediate string
 			bs.WriteByte('>')

pkg/obiformats/fastseq_write_fastq.go

@@ -64,6 +64,9 @@ func FormatFastqBatch(batch obiiter.BioSequenceBatch,
 	first := true
 
 	for _, seq := range batch.Slice() {
+		if len(seq.Id()) == 0 {
+			log.Fatalf("Sequence identifier is empty")
+		}
 		if seq.Len() > 0 {
 			_formatFastq(&bs, seq, formater)
 

pkg/obiformats/json_read.go

@@ -0,0 +1,147 @@
+package obiformats
+
+import (
+	"io"
+	"os"
+	"path"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiiter"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiutils"
+	"github.com/buger/jsonparser"
+	"github.com/goccy/go-json"
+	log "github.com/sirupsen/logrus"
+)
+
+// _parse_json_record parses a single JSON object describing a sequence
+// (as produced by JSONRecord in json_writer.go) into a *obiseq.BioSequence.
+func _parse_json_record(raw []byte, shift byte) *obiseq.BioSequence {
+	sequence := obiseq.NewEmptyBioSequence(0)
+
+	if id, err := jsonparser.GetString(raw, "id"); err == nil {
+		sequence.SetId(id)
+	}
+
+	if seq, err := jsonparser.GetString(raw, "sequence"); err == nil {
+		sequence.SetSequence([]byte(seq))
+	}
+
+	if qual, err := jsonparser.GetString(raw, "qualities"); err == nil {
+		q := []byte(qual)
+		for i := 0; i < len(q); i++ {
+			q[i] -= shift
+		}
+		sequence.SetQualities(q)
+	}
+
+	if annot, dataType, _, err := jsonparser.Get(raw, "annotations"); err == nil && dataType == jsonparser.Object {
+		jsonparser.ObjectEach(annot,
+			func(key []byte, value []byte, valType jsonparser.ValueType, offset int) error {
+				return _parse_json_annotation_field(key, value, valType, sequence)
+			},
+		)
+	}
+
+	return sequence
+}
+
+// _ParseJsonFile streams the top-level JSON array, decoding and pushing one
+// batch of sequences at a time, without ever loading the whole document in
+// memory. Only one raw record at a time is buffered by the decoder.
+func _ParseJsonFile(source string,
+	reader io.Reader,
+	out obiiter.IBioSequence,
+	shift byte,
+	batchSize int) {
+
+	dec := json.NewDecoder(reader)
+
+	if _, err := dec.Token(); err != nil {
+		if err == io.EOF {
+			out.Done()
+			return
+		}
+		log.Fatalf("cannot parse JSON data: %v", err)
+	}
+
+	slice := obiseq.MakeBioSequenceSlice()
+	o := 0
+
+	for dec.More() {
+		var raw json.RawMessage
+
+		if err := dec.Decode(&raw); err != nil {
+			log.Fatalf("cannot parse JSON data: %v", err)
+		}
+
+		sequence := _parse_json_record(raw, shift)
+
+		slice = append(slice, sequence)
+		if len(slice) >= batchSize {
+			out.Push(obiiter.MakeBioSequenceBatch(source, o, slice))
+			o++
+			slice = obiseq.MakeBioSequenceSlice()
+		}
+	}
+
+	if len(slice) > 0 {
+		out.Push(obiiter.MakeBioSequenceBatch(source, o, slice))
+	}
+
+	out.Done()
+}
+
+func ReadJSON(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, error) {
+
+	opt := MakeOptions(options)
+	out := obiiter.MakeIBioSequence()
+
+	out.Add(1)
+	go _ParseJsonFile(opt.Source(),
+		reader,
+		out,
+		obidefault.ReadQualitiesShift(),
+		opt.BatchSize())
+
+	go func() {
+		out.WaitAndClose()
+	}()
+
+	return out, nil
+
+}
+
+func ReadJSONFromFile(filename string, options ...WithOption) (obiiter.IBioSequence, error) {
+
+	options = append(options, OptionsSource(obiutils.RemoveAllExt((path.Base(filename)))))
+	file, err := obiutils.Ropen(filename)
+
+	if err == obiutils.ErrNoContent {
+		log.Infof("file %s is empty", filename)
+		return ReadEmptyFile(options...)
+	}
+
+	if err != nil {
+		return obiiter.NilIBioSequence, err
+	}
+
+	return ReadJSON(file, options...)
+}
+
+func ReadJSONFromStdin(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, error) {
+	options = append(options, OptionsSource(obiutils.RemoveAllExt("stdin")))
+	input, err := obiutils.Buf(os.Stdin)
+
+	if err == obiutils.ErrNoContent {
+		log.Infof("stdin is empty")
+		return ReadEmptyFile(options...)
+	}
+
+	if err != nil {
+		log.Fatalf("open file error: %v", err)
+		return obiiter.NilIBioSequence, err
+	}
+
+	return ReadJSON(input, options...)
+}

pkg/obiformats/ngsfilter_read.go

@@ -631,9 +631,9 @@ func ReadCSVNGSFilter(reader io.Reader) (*obingslibrary.NGSLibrary, error) {
 			return nil, fmt.Errorf("row %d has %d columns, expected %d", len(data), len(fields), len(header))
 		}
 
-		forward_primer := fields[forward_primerColIndex]
+		forward_primer := strings.TrimSpace(fields[forward_primerColIndex])
-		reverse_primer := fields[reverse_primerColIndex]
+		reverse_primer := strings.TrimSpace(fields[reverse_primerColIndex])
-		tags := _parseMainNGSFilterTags(fields[sample_tagColIndex])
+		tags := _parseMainNGSFilterTags(strings.TrimSpace(fields[sample_tagColIndex]))
 
 		marker, _ := ngsfilter.GetMarker(forward_primer, reverse_primer)
 		pcr, ok := marker.GetPCR(tags.Forward, tags.Reverse)
@@ -644,8 +644,8 @@ func ReadCSVNGSFilter(reader io.Reader) (*obingslibrary.NGSLibrary, error) {
 					i, tags.Forward, tags.Reverse, forward_primer, reverse_primer)
 		}
 
-		pcr.Experiment = fields[experimentColIndex]
+		pcr.Experiment = strings.TrimSpace(fields[experimentColIndex])
-		pcr.Sample = fields[sampleColIndex]
+		pcr.Sample = strings.TrimSpace(fields[sampleColIndex])
 
 		if extraColumns != nil {
 			pcr.Annotations = make(obiseq.Annotation)

pkg/obiformats/options.go

@@ -35,6 +35,7 @@ type __options__ struct {
 	csv_auto              bool
 	paired_filename       string
 	source                string
+	filename              string
 	with_feature_table    bool
 	with_pattern          bool
 	with_parent           bool
@@ -216,6 +217,16 @@ func (opt Options) Source() string {
 	return opt.pointer.source
 }
 
+// FileName returns the full path of the file being read, for use in
+// diagnostic messages. It falls back to Source() when no explicit
+// file name has been set (e.g. reading from stdin or a raw reader).
+func (opt Options) FileName() string {
+	if opt.pointer.filename == "" {
+		return opt.pointer.source
+	}
+	return opt.pointer.filename
+}
+
 func (opt Options) WithFeatureTable() bool {
 	return opt.pointer.with_feature_table
 }
@@ -421,6 +432,14 @@ func OptionsSource(source string) WithOption {
 	return f
 }
 
+func OptionsFileName(filename string) WithOption {
+	f := WithOption(func(opt Options) {
+		opt.pointer.filename = filename
+	})
+
+	return f
+}
+
 func OptionsWithProgressBar() WithOption {
 	f := WithOption(func(opt Options) {
 		opt.pointer.with_progress_bar = true

pkg/obiformats/universal_read.go

@@ -145,6 +145,8 @@ func ReadSequencesFromFile(filename string,
 		return ReadGenbank(reader, options...)
 	case "text/csv":
 		return ReadCSV(reader, options...)
+	case "application/json":
+		return ReadJSON(reader, options...)
 	default:
 		log.Fatalf("File %s has guessed format %s which is not yet implemented",
 			filename, mime.String())

pkg/obifp/uint128_test.go

@@ -134,7 +134,7 @@ func TestUint128_QuoRem(t *testing.T) {
 	u := Uint128{w1: 3, w0: 8}
 	v := Uint128{w1: 0, w0: 4}
 	q, r := u.QuoRem(v)
-	assert.Equal(t, Uint128{w1: 0, w0: 2}, q)
+	assert.Equal(t, Uint128{w1: 0, w0: 13835058055282163714}, q)
 	assert.Equal(t, Uint128{w1: 0, w0: 0}, r)
 }
 
@@ -150,7 +150,7 @@ func TestUint128_Div(t *testing.T) {
 	u := Uint128{w1: 3, w0: 8}
 	v := Uint128{w1: 0, w0: 4}
 	q := u.Div(v)
-	assert.Equal(t, Uint128{w1: 0, w0: 2}, q)
+	assert.Equal(t, Uint128{w1: 0, w0: 13835058055282163714}, q)
 }
 
 func TestUint128_Div64(t *testing.T) {
@@ -183,7 +183,7 @@ func TestUint128_Cmp(t *testing.T) {
 func TestUint128_Cmp64(t *testing.T) {
 	u := Uint128{w1: 1, w0: 2}
 	v := uint64(3)
-	assert.Equal(t, -1, u.Cmp64(v))
+	assert.Equal(t, 1, u.Cmp64(v))
 }
 
 func TestUint128_Equals(t *testing.T) {

pkg/obiiter/distribute.go

@@ -57,34 +57,21 @@ func (dist *IDistribute) Classifier() *obiseq.BioSequenceClassifier {
 }
 
 // Distribute organizes the biosequences from the iterator into batches
-// based on the provided classifier and batch sizes. It returns an
+// based on the provided classifier. It returns an IDistribute instance
-// IDistribute instance that manages the distribution of the sequences.
+// that manages the distribution of the sequences.
 //
-// Parameters:
+// Batches are flushed when either BatchSizeMax() sequences or BatchMem()
-//   - class: A pointer to a BioSequenceClassifier used to classify
+// bytes are accumulated per key, mirroring the RebatchBySize strategy.
-//     the biosequences during distribution.
+func (iterator IBioSequence) Distribute(class *obiseq.BioSequenceClassifier) IDistribute {
-//   - sizes: Optional integer values specifying the batch size. If
+	maxCount := obidefault.BatchSizeMax()
-//     no sizes are provided, a default batch size of 5000 is used.
+	maxBytes := obidefault.BatchMem()
-//
-// Returns:
-// An IDistribute instance that contains the outputs of the
-// classified biosequences, a channel for new data notifications,
-// and the classifier used for distribution. The method operates
-// asynchronously, processing the sequences in separate goroutines.
-// It ensures that the outputs are closed and cleaned up once
-// processing is complete.
-func (iterator IBioSequence) Distribute(class *obiseq.BioSequenceClassifier, sizes ...int) IDistribute {
-	batchsize := obidefault.BatchSize()
 
 	outputs := make(map[int]IBioSequence, 100)
 	slices := make(map[int]*obiseq.BioSequenceSlice, 100)
+	bufBytes := make(map[int]int, 100)
 	orders := make(map[int]int, 100)
 	news := make(chan int)
 
-	if len(sizes) > 0 {
-		batchsize = sizes[0]
-	}
-
 	jobDone := sync.WaitGroup{}
 	lock := sync.Mutex{}
 
@@ -115,6 +102,7 @@ func (iterator IBioSequence) Distribute(class *obiseq.BioSequenceClassifier, siz
 					slice = &s
 					slices[key] = slice
 					orders[key] = 0
+					bufBytes[key] = 0
 
 					lock.Lock()
 					outputs[key] = MakeIBioSequence()
@@ -123,14 +111,20 @@ func (iterator IBioSequence) Distribute(class *obiseq.BioSequenceClassifier, siz
 					news <- key
 				}
 
-				*slice = append(*slice, s)
+				sz := s.MemorySize()
-
+				countFull := maxCount > 0 && len(*slice) >= maxCount
-				if len(*slice) == batchsize {
+				memFull := maxBytes > 0 && bufBytes[key]+sz > maxBytes && len(*slice) > 0
+				if countFull || memFull {
 					outputs[key].Push(MakeBioSequenceBatch(source, orders[key], *slice))
 					orders[key]++
 					s := obiseq.MakeBioSequenceSlice()
 					slices[key] = &s
+					slice = &s
+					bufBytes[key] = 0
 				}
+
+				*slice = append(*slice, s)
+				bufBytes[key] += sz
 			}
 		}
 

pkg/obikmer/encodefourmer.go

@@ -47,7 +47,7 @@ func Encode4mer(seq *obiseq.BioSequence, buffer *[]byte) []byte {
 	length := slength - 3
 	rawseq := seq.Sequence()
 
-	if length < 0 {
+	if length <= 0 {
 		return nil
 	}
 

pkg/obikmer/entropy.go

@@ -4,22 +4,21 @@ import "math"
 
 // KmerEntropy computes the entropy of a single encoded k-mer.
 //
-// The algorithm mirrors the lowmask entropy calculation: it decodes the k-mer
+// The algorithm mirrors the Rust obiskbuilder entropy: it decodes the k-mer
 // to a DNA sequence, extracts all sub-words of each size from 1 to levelMax,
 // normalizes them by circular canonical form, counts their frequencies, and
-// computes Shannon entropy normalized by the maximum possible entropy.
+// computes Shannon entropy corrected for class sizes, normalized by the
+// maximum possible entropy over 4^ws raw bins.
 // The returned value is the minimum entropy across all word sizes.
 //
+// Correction for small sequences: the raw entropy H = log(N) - Σ f·log(f)/N
+// under-estimates the true complexity when many raw words collapse to the same
+// canonical form.  Adding Σ f·log(class_size)/N recovers the entropy of the
+// underlying uncollapsed distribution (assuming uniform mixing within each
+// equivalence class).
+//
 // A value close to 0 indicates very low complexity (e.g. "AAAA..."),
 // while a value close to 1 indicates high complexity.
-//
-// Parameters:
-//   - kmer: the encoded k-mer (2 bits per base)
-//   - k: the k-mer size
-//   - levelMax: maximum sub-word size for entropy (typically 6)
-//
-// Returns:
-//   - minimum normalized entropy across all word sizes 1..levelMax
 func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
 	if k < 1 || levelMax < 1 {
 		return 1.0
@@ -35,7 +34,7 @@ func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
 	var seqBuf [32]byte
 	seq := DecodeKmer(kmer, k, seqBuf[:])
 
-	// Pre-compute nLogN lookup (same as lowmask)
+	// Pre-compute nLogN lookup
 	nLogN := make([]float64, k+1)
 	for i := 1; i <= k; i++ {
 		nLogN[i] = float64(i) * math.Log(float64(i))
@@ -51,6 +50,23 @@ func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
 		}
 	}
 
+	// Build ln(class_size) tables: for each canonical form, how many raw
+	// words map to it under circular normalization.
+	classLogSizeTables := make([][]float64, levelMax+1)
+	for ws := 1; ws <= levelMax; ws++ {
+		tableSize := 1 << (ws * 2)
+		classSize := make([]int, tableSize)
+		for code := 0; code < tableSize; code++ {
+			classSize[normTables[ws][code]]++
+		}
+		classLogSizeTables[ws] = make([]float64, tableSize)
+		for j := 0; j < tableSize; j++ {
+			if classSize[j] > 0 {
+				classLogSizeTables[ws][j] = math.Log(float64(classSize[j]))
+			}
+		}
+	}
+
 	minEntropy := math.MaxFloat64
 
 	for ws := 1; ws <= levelMax; ws++ {
@@ -75,23 +91,13 @@ func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
 			table[normWord]++
 		}
 
-		// Compute Shannon entropy
+		// Compute emax over 4^ws raw bins (uncollapsed distribution).
 		floatNwords := float64(nwords)
 		logNwords := math.Log(floatNwords)
-
+		na := tableSize // 4^ws
-		var sumNLogN float64
-		for j := 0; j < tableSize; j++ {
-			n := table[j]
-			if n > 0 {
-				sumNLogN += nLogN[n]
-			}
-		}
-
-		// Compute emax (maximum possible entropy for this word size)
-		na := CanonicalCircularKmerCount(ws)
 		var emax float64
 		if nwords < na {
-			emax = math.Log(float64(nwords))
+			emax = logNwords
 		} else {
 			cov := nwords / na
 			remains := nwords - (na * cov)
@@ -105,7 +111,19 @@ func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
 			continue
 		}
 
-		entropy := (logNwords - sumNLogN/floatNwords) / emax
+		// Accumulate Σ f·log(f) and Σ f·log(class_size) over canonical forms.
+		classLogSize := classLogSizeTables[ws]
+		var sumNLogN, sumClassLogN float64
+		for j := 0; j < tableSize; j++ {
+			n := table[j]
+			if n > 0 {
+				sumNLogN += nLogN[n]
+				sumClassLogN += float64(n) * classLogSize[j]
+			}
+		}
+
+		// Corrected entropy: H_raw ≈ log(N) + (Σf·log(s) - Σf·log(f)) / N
+		entropy := (logNwords + sumClassLogN/floatNwords - sumNLogN/floatNwords) / emax
 		if entropy < 0 {
 			entropy = 0
 		}
@@ -134,6 +152,7 @@ type KmerEntropyFilter struct {
 	threshold           float64
 	nLogN               []float64
 	normTables          [][]int
+	classLogSizeTables  [][]float64
 	emaxValues          []float64
 	logNwords           []float64
 	// Pre-allocated frequency tables reused across Entropy() calls.
@@ -142,11 +161,6 @@ type KmerEntropyFilter struct {
 }
 
 // NewKmerEntropyFilter creates an entropy filter with pre-computed tables.
-//
-// Parameters:
-//   - k: the k-mer size
-//   - levelMax: maximum sub-word size for entropy (typically 6)
-//   - threshold: entropy threshold (k-mers with entropy <= threshold are rejected)
 func NewKmerEntropyFilter(k, levelMax int, threshold float64) *KmerEntropyFilter {
 	if levelMax >= k {
 		levelMax = k - 1
@@ -169,20 +183,38 @@ func NewKmerEntropyFilter(k, levelMax int, threshold float64) *KmerEntropyFilter
 		}
 	}
 
+	// ln(class_size) for each canonical form under circular normalization.
+	classLogSizeTables := make([][]float64, levelMax+1)
+	for ws := 1; ws <= levelMax; ws++ {
+		tableSize := 1 << (ws * 2)
+		classSize := make([]int, tableSize)
+		for code := 0; code < tableSize; code++ {
+			classSize[normTables[ws][code]]++
+		}
+		classLogSizeTables[ws] = make([]float64, tableSize)
+		for j := 0; j < tableSize; j++ {
+			if classSize[j] > 0 {
+				classLogSizeTables[ws][j] = math.Log(float64(classSize[j]))
+			}
+		}
+	}
+
+	// Pre-compute emax and logNwords per word size.
+	// emax uses 4^ws raw bins to match the corrected entropy.
 	emaxValues := make([]float64, levelMax+1)
 	logNwords := make([]float64, levelMax+1)
 	for ws := 1; ws <= levelMax; ws++ {
 		nw := k - ws + 1
-		na := CanonicalCircularKmerCount(ws)
+		na := 1 << (ws * 2) // 4^ws raw bins
+		floatNw := float64(nw)
+		logNwords[ws] = math.Log(floatNw)
 		if nw < na {
-			logNwords[ws] = math.Log(float64(nw))
+			emaxValues[ws] = logNwords[ws]
-			emaxValues[ws] = math.Log(float64(nw))
 		} else {
 			cov := nw / na
 			remains := nw - (na * cov)
-			f1 := float64(cov) / float64(nw)
+			f1 := float64(cov) / floatNw
-			f2 := float64(cov+1) / float64(nw)
+			f2 := float64(cov+1) / floatNw
-			logNwords[ws] = math.Log(float64(nw))
 			emaxValues[ws] = -(float64(na-remains)*f1*math.Log(f1) +
 				float64(remains)*f2*math.Log(f2))
 		}
@@ -200,6 +232,7 @@ func NewKmerEntropyFilter(k, levelMax int, threshold float64) *KmerEntropyFilter
 		threshold:          threshold,
 		nLogN:              nLogN,
 		normTables:         normTables,
+		classLogSizeTables: classLogSizeTables,
 		emaxValues:         emaxValues,
 		logNwords:          logNwords,
 		freqTables:         freqTables,
@@ -236,7 +269,7 @@ func (ef *KmerEntropyFilter) Entropy(kmer uint64) float64 {
 		// Count circular-canonical sub-word frequencies
 		tableSize := 1 << (ws * 2)
 		table := ef.freqTables[ws]
-		clear(table) // reset to zero
+		clear(table)
 		mask := (1 << (ws * 2)) - 1
 		normTable := ef.normTables[ws]
 
@@ -251,19 +284,21 @@ func (ef *KmerEntropyFilter) Entropy(kmer uint64) float64 {
 			table[normWord]++
 		}
 
-		// Compute Shannon entropy
 		floatNwords := float64(nwords)
 		logNwords := ef.logNwords[ws]
+		classLogSize := ef.classLogSizeTables[ws]
 
-		var sumNLogN float64
+		var sumNLogN, sumClassLogN float64
 		for j := 0; j < tableSize; j++ {
 			n := table[j]
 			if n > 0 {
 				sumNLogN += ef.nLogN[n]
+				sumClassLogN += float64(n) * classLogSize[j]
 			}
 		}
 
-		entropy := (logNwords - sumNLogN/floatNwords) / emax
+		// Corrected entropy: H_raw ≈ log(N) + (Σf·log(s) - Σf·log(f)) / N
+		entropy := (logNwords + sumClassLogN/floatNwords - sumNLogN/floatNwords) / emax
 		if entropy < 0 {
 			entropy = 0
 		}

pkg/obilua/lua.go

@@ -91,7 +91,7 @@ func LuaWorker(proto *lua.FunctionProto) obiseq.SeqWorker {
 	err := interpreter.PCall(0, lua.MultRet, nil)
 
 	if err != nil {
-		log.Fatalf("Error in executing the lua script")
+		log.Fatalf("Error in executing the lua script: %v", err)
 	}
 
 	result := interpreter.GetGlobal("worker")
@@ -141,6 +141,69 @@ func LuaWorker(proto *lua.FunctionProto) obiseq.SeqWorker {
 	return nil
 }
 
+// LuaSliceWorker creates a SeqSliceWorker that calls the Lua function
+// named "slice_worker". Unlike LuaWorker, the entire batch (BioSequenceSlice)
+// is passed to the Lua function at once, enabling batch-level processing
+// (e.g. a single HTTP request per batch instead of one per sequence).
+//
+// The Lua function signature:
+//
+//	function slice_worker(slice)   -- receives a BioSequenceSlice
+//	    -- process the batch
+//	    return slice               -- returns a BioSequenceSlice (or nil)
+//	end
+func LuaSliceWorker(proto *lua.FunctionProto) obiseq.SeqSliceWorker {
+	interpreter := NewInterpreter()
+	lfunc := interpreter.NewFunctionFromProto(proto)
+	interpreter.Push(lfunc)
+	err := interpreter.PCall(0, lua.MultRet, nil)
+
+	if err != nil {
+		log.Fatalf("Error in executing the lua script: %v", err)
+	}
+
+	result := interpreter.GetGlobal("slice_worker")
+
+	if lua_worker, ok := result.(*lua.LFunction); ok {
+		f := func(slice obiseq.BioSequenceSlice) (obiseq.BioSequenceSlice, error) {
+			if err := interpreter.CallByParam(lua.P{
+				Fn:      lua_worker,
+				NRet:    1,
+				Protect: true,
+			}, obiseqslice2Lua(interpreter, &slice)); err != nil {
+				log.Fatal(err)
+			}
+
+			lreponse := interpreter.Get(-1)
+			defer interpreter.Pop(1)
+
+			if reponse, ok := lreponse.(*lua.LUserData); ok {
+				s := reponse.Value
+				switch val := s.(type) {
+				case *obiseq.BioSequenceSlice:
+					return *val, nil
+				case *obiseq.BioSequence:
+					return obiseq.BioSequenceSlice{val}, nil
+				default:
+					r := reflect.TypeOf(val)
+					return nil, fmt.Errorf("slice_worker function doesn't return the correct type %s", r)
+				}
+			}
+
+			if _, ok = lreponse.(*lua.LNilType); ok {
+				return nil, nil
+			}
+
+			return nil, fmt.Errorf("slice_worker function doesn't return the correct type %T", lreponse)
+		}
+
+		return f
+	}
+
+	log.Fatalf("The slice_worker object is not a function")
+	return nil
+}
+
 // LuaProcessor processes a Lua script on a sequence iterator and returns a new iterator.
 //
 // Parameters:
@@ -173,7 +236,7 @@ func LuaProcessor(iterator obiiter.IBioSequence, name, program string, breakOnEr
 	err = interpreter.PCall(0, lua.MultRet, nil)
 
 	if err != nil {
-		log.Fatalf("Error in executing the lua script")
+		log.Fatalf("Error in executing the lua script: %v", err)
 	}
 
 	result := interpreter.GetGlobal("begin")
@@ -198,7 +261,7 @@ func LuaProcessor(iterator obiiter.IBioSequence, name, program string, breakOnEr
 		err = interpreter.PCall(0, lua.MultRet, nil)
 
 		if err != nil {
-			log.Fatalf("Error in executing the lua script")
+			log.Fatalf("Error in executing the lua script: %v", err)
 		}
 
 		result := interpreter.GetGlobal("finish")
@@ -216,11 +279,27 @@ func LuaProcessor(iterator obiiter.IBioSequence, name, program string, breakOnEr
 
 	}()
 
-	ff := func(iterator obiiter.IBioSequence) {
+	// Detect whether the script defines slice_worker (batch-level) or worker (per-sequence).
-		w := LuaWorker(proto)
+	hasSliceWorker := func() bool {
-		sw := obiseq.SeqToSliceWorker(w, false)
+		interpreter := NewInterpreter()
+		lfunc := interpreter.NewFunctionFromProto(proto)
+		interpreter.Push(lfunc)
+		if err := interpreter.PCall(0, lua.MultRet, nil); err != nil {
+			return false
+		}
+		result := interpreter.GetGlobal("slice_worker")
+		interpreter.Close()
+		_, ok := result.(*lua.LFunction)
+		return ok
+	}()
 
-		// iterator = iterator.SortBatches()
+	ff := func(iterator obiiter.IBioSequence) {
+		var sw obiseq.SeqSliceWorker
+		if hasSliceWorker {
+			sw = LuaSliceWorker(proto)
+		} else {
+			sw = obiseq.SeqToSliceWorker(LuaWorker(proto), false)
+		}
 
 		for iterator.Next() {
 			seqs := iterator.Get()
@@ -235,6 +314,10 @@ func LuaProcessor(iterator obiiter.IBioSequence, name, program string, breakOnEr
 				}
 			}
 
+			if ns == nil {
+				ns = obiseq.BioSequenceSlice{}
+			}
+
 			newIter.Push(obiiter.MakeBioSequenceBatch(seqs.Source(), seqs.Order(), ns))
 		}
 

pkg/obilua/lua_push_interface.go

@@ -17,15 +17,7 @@ import (
 // No return values. This function operates directly on the Lua state stack.
 func pushInterfaceToLua(L *lua.LState, val interface{}) {
 	switch v := val.(type) {
-	case string:
+	// Typed slices and maps from internal OBITools code — not produced by json.Unmarshal
-		L.Push(lua.LString(v))
-	case bool:
-		L.Push(lua.LBool(v))
-	case int:
-		L.Push(lua.LNumber(v))
-	case float64:
-		L.Push(lua.LNumber(v))
-	// Add other cases as needed for different types
 	case map[string]int:
 		pushMapStringIntToLua(L, v)
 	case map[string]string:
@@ -34,8 +26,6 @@ func pushInterfaceToLua(L *lua.LState, val interface{}) {
 		pushMapStringBoolToLua(L, v)
 	case map[string]float64:
 		pushMapStringFloat64ToLua(L, v)
-	case map[string]interface{}:
-		pushMapStringInterfaceToLua(L, v)
 	case []string:
 		pushSliceStringToLua(L, v)
 	case []int:
@@ -46,63 +36,63 @@ func pushInterfaceToLua(L *lua.LState, val interface{}) {
 		pushSliceNumericToLua(L, v)
 	case []bool:
 		pushSliceBoolToLua(L, v)
-	case []interface{}:
-		pushSliceInterfaceToLua(L, v)
-	case nil:
-		L.Push(lua.LNil)
 	case *sync.Mutex:
 		pushMutexToLua(L, v)
 	default:
-		log.Fatalf("Cannot deal with value (%T) : %v", val, val)
+		// Handles nil, bool, int, float64, string, map[string]interface{},
+		// []interface{} — all recursively via lvalueFromInterface.
+		L.Push(lvalueFromInterface(L, v))
 	}
 }
 
 func pushMapStringInterfaceToLua(L *lua.LState, m map[string]interface{}) {
-	// Create a new Lua table
 	luaTable := L.NewTable()
-	// Iterate over the Go map and set the key-value pairs in the Lua table
 	for key, value := range m {
-		switch v := value.(type) {
+		L.SetField(luaTable, key, lvalueFromInterface(L, value))
-		case int:
-			luaTable.RawSetString(key, lua.LNumber(v))
-		case float64:
-			luaTable.RawSetString(key, lua.LNumber(v))
-		case bool:
-			luaTable.RawSetString(key, lua.LBool(v))
-		case string:
-			luaTable.RawSetString(key, lua.LString(v))
-		default:
-			log.Fatalf("Doesn't deal with map containing value %v of type %T", v, v)
 	}
-	}
-
-	// Push the Lua table onto the stack
 	L.Push(luaTable)
 }
 
 func pushSliceInterfaceToLua(L *lua.LState, s []interface{}) {
-	// Create a new Lua table
 	luaTable := L.NewTable()
-	// Iterate over the Go map and set the key-value pairs in the Lua table
 	for _, value := range s {
-		switch v := value.(type) {
+		luaTable.Append(lvalueFromInterface(L, value))
-		case int:
-			luaTable.Append(lua.LNumber(v))
-		case float64:
-			luaTable.Append(lua.LNumber(v))
-		case bool:
-			luaTable.Append(lua.LBool(v))
-		case string:
-			luaTable.Append(lua.LString(v))
-		default:
-			log.Fatalf("Doesn't deal with slice containing value %v of type %T", v, v)
 	}
-	}
-
-	// Push the Lua table onto the stack
 	L.Push(luaTable)
 }
 
+// lvalueFromInterface converts a Go interface{} value (as produced by json.Unmarshal)
+// to the corresponding lua.LValue, handling nested maps and slices recursively.
+func lvalueFromInterface(L *lua.LState, value interface{}) lua.LValue {
+	switch v := value.(type) {
+	case nil:
+		return lua.LNil
+	case bool:
+		return lua.LBool(v)
+	case int:
+		return lua.LNumber(v)
+	case float64:
+		return lua.LNumber(v)
+	case string:
+		return lua.LString(v)
+	case map[string]interface{}:
+		t := L.NewTable()
+		for key, val := range v {
+			L.SetField(t, key, lvalueFromInterface(L, val))
+		}
+		return t
+	case []interface{}:
+		t := L.NewTable()
+		for _, val := range v {
+			t.Append(lvalueFromInterface(L, val))
+		}
+		return t
+	default:
+		log.Fatalf("lvalueFromInterface: unsupported type %T: %v", v, v)
+		return lua.LNil
+	}
+}
+
 // pushMapStringIntToLua creates a new Lua table and iterates over the Go map to set key-value pairs in the Lua table. It then pushes the Lua table onto the stack.
 //
 // L *lua.LState - the Lua state

pkg/obilua/lua_table.go

@@ -28,6 +28,8 @@ func Table2Interface(interpreter *lua.LState, table *lua.LTable) interface{} {
 				val[i-1] = float64(v.(lua.LNumber))
 			case lua.LTString:
 				val[i-1] = string(v.(lua.LString))
+			case lua.LTTable:
+				val[i-1] = Table2Interface(interpreter, v.(*lua.LTable))
 			}
 		}
 		return val
@@ -45,6 +47,8 @@ func Table2Interface(interpreter *lua.LState, table *lua.LTable) interface{} {
 					val[string(ks)] = float64(v.(lua.LNumber))
 				case lua.LTString:
 					val[string(ks)] = string(v.(lua.LString))
+				case lua.LTTable:
+					val[string(ks)] = Table2Interface(interpreter, v.(*lua.LTable))
 				}
 			}
 		})

pkg/obilua/luahttp.go

@@ -0,0 +1,128 @@
+package obilua
+
+import (
+	"context"
+	"io"
+	"net/http"
+	"strings"
+	"sync"
+	"time"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault"
+	lua "github.com/yuin/gopher-lua"
+)
+
+const httpClientTimeout = 300 * time.Second
+
+var (
+	_httpClient     *http.Client
+	_httpClientOnce sync.Once
+
+	// _httpSemaphore limits the number of concurrent HTTP requests.
+	// Initialised lazily alongside the client.
+	_httpSemaphore chan struct{}
+)
+
+func getHTTPClient() *http.Client {
+	_httpClientOnce.Do(func() {
+		conns := 2 * obidefault.ParallelWorkers()
+		_httpClient = &http.Client{
+			Transport: &http.Transport{
+				MaxIdleConnsPerHost: conns,
+				MaxConnsPerHost:     conns,
+				IdleConnTimeout:     90 * time.Second,
+			},
+			Timeout: httpClientTimeout,
+		}
+		_httpSemaphore = make(chan struct{}, obidefault.ParallelWorkers())
+	})
+	return _httpClient
+}
+
+// RegisterHTTP registers the http module in the Lua state as a global,
+// consistent with obicontext and BioSequence.
+//
+// Exposes:
+//
+//	http.post(url, body [, timeout_ms]) → response string  (on success)
+//	http.post(url, body [, timeout_ms]) → nil, err string  (on error)
+//	http.set_concurrency(n)             → set max simultaneous requests
+func RegisterHTTP(luaState *lua.LState) {
+	table := luaState.NewTable()
+	luaState.SetField(table, "post", luaState.NewFunction(luaHTTPPost))
+	luaState.SetField(table, "set_concurrency", luaState.NewFunction(luaHTTPSetConcurrency))
+	luaState.SetGlobal("http", table)
+}
+
+// luaHTTPPost implements http.post(url, body [, timeout_ms]) for Lua.
+//
+// The optional third argument overrides the default timeout (in milliseconds).
+// Concurrent requests are throttled through _httpSemaphore so that a
+// single-threaded backend server is not overwhelmed by K parallel workers.
+//
+// Lua signature:
+//
+//	local response          = http.post(url, body)
+//	local response          = http.post(url, body, 5000)   -- 5 s timeout
+//	local response, err     = http.post(url, body)
+func luaHTTPPost(L *lua.LState) int {
+	url := L.CheckString(1)
+	body := L.CheckString(2)
+
+	client := getHTTPClient()
+
+	timeout := httpClientTimeout
+	if L.GetTop() >= 3 {
+		ms := L.CheckInt(3)
+		timeout = time.Duration(ms) * time.Millisecond
+	}
+
+	// Acquire semaphore slot — blocks until a slot is free.
+	_httpSemaphore <- struct{}{}
+	defer func() { <-_httpSemaphore }()
+
+	ctx, cancel := context.WithTimeout(context.Background(), timeout)
+	defer cancel()
+
+	req, err := http.NewRequestWithContext(ctx, http.MethodPost, url, strings.NewReader(body))
+	if err != nil {
+		L.Push(lua.LNil)
+		L.Push(lua.LString(err.Error()))
+		return 2
+	}
+	req.Header.Set("Content-Type", "application/json")
+
+	resp, err := client.Do(req)
+	if err != nil {
+		L.Push(lua.LNil)
+		L.Push(lua.LString(err.Error()))
+		return 2
+	}
+	defer resp.Body.Close()
+
+	respBytes, err := io.ReadAll(resp.Body)
+	if err != nil {
+		L.Push(lua.LNil)
+		L.Push(lua.LString(err.Error()))
+		return 2
+	}
+
+	L.Push(lua.LString(respBytes))
+	return 1
+}
+
+// luaHTTPSetConcurrency replaces the semaphore with a new one of size n.
+// Must be called before the first http.post (e.g. in begin()).
+//
+// Lua signature:
+//
+//	http.set_concurrency(1)   -- serialise all HTTP requests
+func luaHTTPSetConcurrency(L *lua.LState) int {
+	n := L.CheckInt(1)
+	if n < 1 {
+		n = 1
+	}
+	getHTTPClient() // ensure singleton is initialised
+	_httpSemaphore = make(chan struct{}, n)
+	return 0
+}

pkg/obilua/luajson.go

@@ -0,0 +1,71 @@
+package obilua
+
+import (
+	"encoding/json"
+
+	lua "github.com/yuin/gopher-lua"
+)
+
+// RegisterJSON registers the json module in the Lua state as a global,
+// consistent with obicontext, BioSequence, and http.
+//
+// Exposes:
+//
+//	json.encode(data)   → string         (on success)
+//	json.encode(data)   → nil, err       (on error)
+//	json.decode(string) → value          (on success)
+//	json.decode(string) → nil, err       (on error)
+func RegisterJSON(luaState *lua.LState) {
+	table := luaState.NewTable()
+	luaState.SetField(table, "encode", luaState.NewFunction(luaJSONEncode))
+	luaState.SetField(table, "decode", luaState.NewFunction(luaJSONDecode))
+	luaState.SetGlobal("json", table)
+}
+
+// luaJSONEncode implements json.encode(data) for Lua.
+func luaJSONEncode(L *lua.LState) int {
+	val := L.CheckAny(1)
+
+	var goVal interface{}
+	switch v := val.(type) {
+	case *lua.LTable:
+		goVal = Table2Interface(L, v)
+	case lua.LString:
+		goVal = string(v)
+	case lua.LNumber:
+		goVal = float64(v)
+	case lua.LBool:
+		goVal = bool(v)
+	case *lua.LNilType:
+		goVal = nil
+	default:
+		L.Push(lua.LNil)
+		L.Push(lua.LString("json.encode: unsupported type"))
+		return 2
+	}
+
+	b, err := json.Marshal(goVal)
+	if err != nil {
+		L.Push(lua.LNil)
+		L.Push(lua.LString(err.Error()))
+		return 2
+	}
+
+	L.Push(lua.LString(b))
+	return 1
+}
+
+// luaJSONDecode implements json.decode(string) for Lua.
+func luaJSONDecode(L *lua.LState) int {
+	s := L.CheckString(1)
+
+	var goVal interface{}
+	if err := json.Unmarshal([]byte(s), &goVal); err != nil {
+		L.Push(lua.LNil)
+		L.Push(lua.LString(err.Error()))
+		return 2
+	}
+
+	pushInterfaceToLua(L, goVal)
+	return 1
+}

pkg/obilua/luajson_test.go

@@ -0,0 +1,184 @@
+package obilua
+
+import (
+	"testing"
+
+	lua "github.com/yuin/gopher-lua"
+)
+
+// runLua executes a Lua snippet inside a fresh interpreter and returns the
+// LState so the caller can inspect the stack.
+func runLua(t *testing.T, script string) *lua.LState {
+	t.Helper()
+	L := NewInterpreter()
+	if err := L.DoString(script); err != nil {
+		t.Fatalf("Lua error: %v", err)
+	}
+	return L
+}
+
+// TestJSONEncodeScalar verifies that simple scalars are encoded correctly.
+func TestJSONEncodeScalar(t *testing.T) {
+	cases := []struct {
+		script   string
+		expected string
+	}{
+		{`result = json.encode("hello")`, `"hello"`},
+		{`result = json.encode(42)`, `42`},
+		{`result = json.encode(true)`, `true`},
+	}
+
+	for _, tc := range cases {
+		L := runLua(t, tc.script)
+		got := string(L.GetGlobal("result").(lua.LString))
+		if got != tc.expected {
+			t.Errorf("encode(%s): got %q, want %q", tc.script, got, tc.expected)
+		}
+		L.Close()
+	}
+}
+
+// TestJSONEncodeTable verifies that a Lua table (array and map) encodes to JSON.
+func TestJSONEncodeTable(t *testing.T) {
+	L := runLua(t, `result = json.encode({a = 1, b = "x"})`)
+	got := string(L.GetGlobal("result").(lua.LString))
+	// json.Marshal produces deterministic output for maps in Go 1.12+... actually not.
+	// Just check it round-trips via decode instead.
+	L.Close()
+	if got == "" {
+		t.Fatal("encode returned empty string")
+	}
+}
+
+// TestJSONDecodeScalar verifies that JSON scalars decode to the right Lua types.
+func TestJSONDecodeScalar(t *testing.T) {
+	L := runLua(t, `
+		s = json.decode('"hello"')
+		n = json.decode('3.14')
+		b = json.decode('true')
+	`)
+	if s, ok := L.GetGlobal("s").(lua.LString); !ok || string(s) != "hello" {
+		t.Errorf("decode string: got %v", L.GetGlobal("s"))
+	}
+	if n, ok := L.GetGlobal("n").(lua.LNumber); !ok || float64(n) != 3.14 {
+		t.Errorf("decode number: got %v", L.GetGlobal("n"))
+	}
+	if b, ok := L.GetGlobal("b").(lua.LBool); !ok || !bool(b) {
+		t.Errorf("decode bool: got %v", L.GetGlobal("b"))
+	}
+	L.Close()
+}
+
+// TestJSONRoundTripFlat verifies a flat table survives encode → decode.
+func TestJSONRoundTripFlat(t *testing.T) {
+	L := runLua(t, `
+		original = {name = "Homo_sapiens", score = 1.0, valid = true}
+		encoded  = json.encode(original)
+		decoded  = json.decode(encoded)
+	`)
+	decoded, ok := L.GetGlobal("decoded").(*lua.LTable)
+	if !ok {
+		t.Fatal("decoded is not a table")
+	}
+	if v := decoded.RawGetString("name"); string(v.(lua.LString)) != "Homo_sapiens" {
+		t.Errorf("name: got %v", v)
+	}
+	if v := decoded.RawGetString("score"); float64(v.(lua.LNumber)) != 1.0 {
+		t.Errorf("score: got %v", v)
+	}
+	if v := decoded.RawGetString("valid"); !bool(v.(lua.LBool)) {
+		t.Errorf("valid: got %v", v)
+	}
+	L.Close()
+}
+
+// TestJSONRoundTripNested verifies a 3-level nested structure (kmindex response)
+// survives encode → decode with correct values at every level.
+func TestJSONRoundTripNested(t *testing.T) {
+	L := NewInterpreter()
+
+	// Inject the JSON string as a Lua global to avoid quoting issues.
+	L.SetGlobal("kmindex_json", lua.LString(
+		`{"Human":{"query_001":{"Homo_sapiens--GCF_000001405_40":1.0}}}`,
+	))
+
+	if err := L.DoString(`
+		data      = json.decode(kmindex_json)
+		reencoded = json.encode(data)
+		data2     = json.decode(reencoded)
+	`); err != nil {
+		t.Fatalf("Lua error: %v", err)
+	}
+
+	// Navigate data["Human"]["query_001"]["Homo_sapiens--GCF_000001405_40"]
+	data, ok := L.GetGlobal("data").(*lua.LTable)
+	if !ok {
+		t.Fatal("data is not a table")
+	}
+	human, ok := data.RawGetString("Human").(*lua.LTable)
+	if !ok {
+		t.Fatal("data.Human is not a table")
+	}
+	query, ok := human.RawGetString("query_001").(*lua.LTable)
+	if !ok {
+		t.Fatal("data.Human.query_001 is not a table")
+	}
+	score, ok := query.RawGetString("Homo_sapiens--GCF_000001405_40").(lua.LNumber)
+	if !ok || float64(score) != 1.0 {
+		t.Errorf("score: got %v, want 1.0", query.RawGetString("Homo_sapiens--GCF_000001405_40"))
+	}
+
+	// Same check on the re-encoded+decoded version
+	data2, ok := L.GetGlobal("data2").(*lua.LTable)
+	if !ok {
+		t.Fatal("data2 is not a table")
+	}
+	score2 := data2.RawGetString("Human").(*lua.LTable).
+		RawGetString("query_001").(*lua.LTable).
+		RawGetString("Homo_sapiens--GCF_000001405_40").(lua.LNumber)
+	if float64(score2) != 1.0 {
+		t.Errorf("data2 score: got %v, want 1.0", score2)
+	}
+	L.Close()
+}
+
+// TestJSONDecodeArray verifies that a JSON array decodes to a Lua array table.
+func TestJSONDecodeArray(t *testing.T) {
+	L := runLua(t, `arr = json.decode('[1, 2, 3]')`)
+	arr, ok := L.GetGlobal("arr").(*lua.LTable)
+	if !ok {
+		t.Fatal("arr is not a table")
+	}
+	for i, expected := range []float64{1, 2, 3} {
+		v, ok := arr.RawGetInt(i + 1).(lua.LNumber)
+		if !ok || float64(v) != expected {
+			t.Errorf("arr[%d]: got %v, want %v", i+1, arr.RawGetInt(i+1), expected)
+		}
+	}
+	L.Close()
+}
+
+// TestJSONEncodeError verifies that json.encode on an unsupported type returns nil + error.
+func TestJSONEncodeError(t *testing.T) {
+	L := runLua(t, `
+		local result, err = json.encode(nil)
+	`)
+	// nil encodes to JSON "null" — not an error
+	L.Close()
+}
+
+// TestJSONDecodeError verifies that malformed JSON returns nil + error string.
+func TestJSONDecodeError(t *testing.T) {
+	L := runLua(t, `
+		local result, err = json.decode("not valid json")
+		decode_ok     = (result == nil)
+		decode_has_err = (err ~= nil)
+	`)
+	if L.GetGlobal("decode_ok") != lua.LTrue {
+		t.Error("expected nil result on decode error")
+	}
+	if L.GetGlobal("decode_has_err") != lua.LTrue {
+		t.Error("expected error string on decode error")
+	}
+	L.Close()
+}

pkg/obilua/obilib.go

@@ -5,4 +5,6 @@ import lua "github.com/yuin/gopher-lua"
 func RegisterObilib(luaState *lua.LState) {
 	RegisterObiSeq(luaState)
 	RegisterObiTaxonomy(luaState)
+	RegisterHTTP(luaState)
+	RegisterJSON(luaState)
 }

pkg/obilua/obiseqslice.go

@@ -31,7 +31,8 @@ func obiseqslice2Lua(interpreter *lua.LState,
 }
 
 func newObiSeqSlice(luaState *lua.LState) int {
-	seqslice := obiseq.NewBioSequenceSlice()
+	capacity := luaState.OptInt(1, 0)
+	seqslice := obiseq.NewBioSequenceSlice(capacity)
 	luaState.Push(obiseqslice2Lua(luaState, seqslice))
 	return 1
 }

pkg/obingslibrary/multimatch.go

@@ -777,7 +777,7 @@ func (library *NGSLibrary) ExtractMultiBarcodeSliceWorker(options ...WithOption)
 		library.SetAllowsIndels(true)
 	}
 
-	if opt.AllowedMismatches() > 0 {
+	if opt.AllowedMismatchesIsSet() {
 		library.SetAllowedMismatches(opt.AllowedMismatches())
 	}
 

pkg/obingslibrary/worker.go

@@ -9,6 +9,7 @@ type _Options struct {
 	discardErrors      bool
 	unidentified       string
 	allowedMismatch    int
+	allowedMismatchSet bool
 	allowsIndel        bool
 	withProgressBar    bool
 	parallelWorkers    int
@@ -52,6 +53,7 @@ func OptionWithProgressBar(yes bool) WithOption {
 func OptionAllowedMismatches(count int) WithOption {
 	f := WithOption(func(opt Options) {
 		opt.pointer.allowedMismatch = count
+		opt.pointer.allowedMismatchSet = true
 	})
 
 	return f
@@ -97,6 +99,12 @@ func (options Options) AllowedMismatches() int {
 	return options.pointer.allowedMismatch
 }
 
+// AllowedMismatchesIsSet returns true if OptionAllowedMismatches
+// was explicitly applied to these options.
+func (options Options) AllowedMismatchesIsSet() bool {
+	return options.pointer.allowedMismatchSet
+}
+
 func (options Options) AllowsIndels() bool {
 	return options.pointer.allowsIndel
 }

pkg/obioptions/version.go

@@ -3,7 +3,7 @@ package obioptions
 // Version is automatically updated by the Makefile from version.txt
 // The patch number (third digit) is incremented on each push to the repository
 
-var _Version = "Release 4.4.24"
+var _Version = "Release 4.5.0"
 
 // Version returns the version of the obitools package.
 //

pkg/obiseq/attributes.go

@@ -364,6 +364,24 @@ func (s *BioSequence) GetIntSlice(key string) ([]int, bool) {
 	return val, ok
 }
 
+func (s *BioSequence) GetMapOfIntSlice(key string) (map[string][]int, bool) {
+	v, ok := s.GetAttribute(key)
+	if !ok {
+		return nil, false
+	}
+	val, err := obiutils.InterfaceToMapOfIntSlice(v)
+	return val, err == nil
+}
+
+func (s *BioSequence) GetMapOfStringSlice(key string) (map[string][]string, bool) {
+	v, ok := s.GetAttribute(key)
+	if !ok {
+		return nil, false
+	}
+	val, err := obiutils.InterfaceToMapOfStringSlice(v)
+	return val, err == nil
+}
+
 // Count returns the value of the "count" attribute of the BioSequence.
 //
 // The count of a sequence is the number of times it has been observed in the dataset.

pkg/obiseq/biosequence.go

@@ -499,6 +499,9 @@ func (s *BioSequence) SetQualities(qualities Quality) {
 	if s.qualities != nil {
 		RecycleSlice(&s.qualities)
 	}
+	if len(qualities) > 0 && len(qualities) != len(s.sequence) {
+		log.Panicf("[BioSequence.SetQualities] Sequence %s has a length of %d and qualities a length of %d", s.id, len(s.sequence), len(qualities))
+	}
 	s.qualities = CopySlice(qualities)
 }
 
@@ -508,6 +511,9 @@ func (s *BioSequence) TakeQualities(qualities Quality) {
 	if s.qualities != nil {
 		RecycleSlice(&s.qualities)
 	}
+	if len(qualities) > 0 && len(qualities) != len(s.sequence) {
+		log.Panicf("[BioSequence.TakeQualities] Sequence %s has a length of %d and qualities a length of %d", s.id, len(s.sequence), len(qualities))
+	}
 	s.qualities = qualities
 }
 

pkg/obiseq/biosequence_test.go

@@ -103,7 +103,7 @@ func TestNewBioSequence(t *testing.T) {
 // Return type: None.
 func TestNewBioSequenceWithQualities(t *testing.T) {
 	id := "123"
-	sequence := []byte("ATGC")
+	sequence := []byte("atgc")
 	definition := "DNA sequence"
 	qualities := []byte("1234")
 

pkg/obiseq/language.go

@@ -141,6 +141,33 @@ var OBILang = gval.NewLanguage(
 	gval.Function("max", func(args ...interface{}) (interface{}, error) {
 		return obiutils.Max(args[0])
 	}),
+	gval.Function("whichmax", func(args ...interface{}) (interface{}, error) {
+		result, err := obiutils.WhichMax(args[0])
+		if idx, ok := result.(int); ok {
+			return float64(idx), nil
+		}
+		return result, err
+	}),
+	gval.Function("whichmin", func(args ...interface{}) (interface{}, error) {
+		result, err := obiutils.WhichMin(args[0])
+		if idx, ok := result.(int); ok {
+			return float64(idx), nil
+		}
+		return result, err
+	}),
+
+	gval.Function("filtermin", func(args ...interface{}) (interface{}, error) {
+		return obiutils.FilterMin(args[0], args[1])
+	}),
+
+	gval.Function("filtermax", func(args ...interface{}) (interface{}, error) {
+		return obiutils.FilterMax(args[0], args[1])
+	}),
+
+	gval.Function("saturatingsub", func(args ...interface{}) (interface{}, error) {
+		return obiutils.SaturatingSub(args[0], args[1])
+	}),
+
 	gval.Function("contains", func(args ...interface{}) (interface{}, error) {
 		if obiutils.IsAMap(args[0]) {
 			val := reflect.ValueOf(args[0]).MapIndex(reflect.ValueOf(args[1]))

pkg/obiseq/revcomp.go

@@ -118,6 +118,9 @@ func (sequence *BioSequence) _revcmpMutation() *BioSequence {
  */
 func ReverseComplementWorker(inplace bool) SeqWorker {
 	f := func(input *BioSequence) (BioSequenceSlice, error) {
+		if input.IsPaired() {
+			input.PairedWith().ReverseComplement(inplace)
+		}
 		return BioSequenceSlice{input.ReverseComplement(inplace)}, nil
 	}
 

pkg/obiseq/subseq.go

@@ -48,7 +48,16 @@ func (sequence *BioSequence) Subsequence(from, to int, circular bool) (*BioSeque
 		newSeq.sequence = CopySlice(sequence.Sequence()[from:to])
 
 		if sequence.HasQualities() {
-			newSeq.qualities = CopySlice(sequence.Qualities()[from:to])
+			qual := sequence.Qualities()
+			if len(qual) != sequence.Len() {
+				log.Panicf(
+					"[BioSequence.Subsequence] Sequence %s has a length of %d and qualities a length of %d",
+					sequence.Id(),
+					sequence.Len(),
+					len(qual),
+				)
+			}
+			newSeq.qualities = CopySlice(qual[from:to])
 		}
 
 		newSeq.id = fmt.Sprintf("%s_sub[%d..%d]", sequence.Id(), from+1, to)
@@ -58,7 +67,16 @@ func (sequence *BioSequence) Subsequence(from, to int, circular bool) (*BioSeque
 		newSeq.Write(sequence.Sequence()[0:to])
 
 		if sequence.HasQualities() {
-			newSeq.WriteQualities(sequence.Qualities()[0:to])
+			qual := sequence.Qualities()
+			if len(qual) != sequence.Len() {
+				log.Panicf(
+					"[BioSequence.Subsequence] Sequence %s has a length of %d and qualities a length of %d",
+					sequence.Id(),
+					sequence.Len(),
+					len(qual),
+				)
+			}
+			newSeq.WriteQualities(qual[0:to])
 		}
 
 	}

pkg/obiseq/taxonomy_methods.go

@@ -70,6 +70,12 @@ func (s *BioSequence) SetTaxid(taxid string, rank ...string) {
 				}
 			}
 
+		} else if obidefault.UseRawTaxids() {
+			// Without a loaded taxonomy, extract the bare ID from full-format strings
+			// like "code:12345 [Name]@rank" so that --raw-taxid is honoured everywhere.
+			if _, rawID, _, _, parseErr := obitax.ParseTaxonString(taxid); parseErr == nil {
+				taxid = rawID
+			}
 		}
 	}
 
@@ -177,7 +183,7 @@ func (sequence *BioSequence) SetPath(taxonomy *obitax.Taxonomy) []string {
 	lpath := path.Len() - 1
 
 	for i := lpath; i >= 0; i-- {
-		spath[lpath-i] = path.Get(i).String(taxonomy.Code())
+		spath[lpath-i] = path.Get(i).FullString(taxonomy.Code())
 	}
 
 	sequence.SetAttribute("taxonomic_path", spath)

pkg/obiseq/worker.go

@@ -104,11 +104,11 @@ func SeqToSliceWorker(worker SeqWorker,
 			for _, s := range input {
 				r, err := worker(s)
 				if err == nil {
-					for _, rs := range r {
+					if i+len(r) > cap(output) {
-						if i == len(output) {
+						output = slices.Grow(output[:i], len(r))
-							output = slices.Grow(output, cap(output))
 						output = output[:cap(output)]
 					}
+					for _, rs := range r {
 						output[i] = rs
 						i++
 					}

pkg/obitax/taxonnode.go

@@ -29,6 +29,24 @@ type TaxNode struct {
 	alternatenames *map[*string]*string
 }
 
+// FullString returns the full string representation of the TaxNode in the form
+// "taxonomyCode:id [scientificName]@rank", regardless of the UseRawTaxids setting.
+// This is used internally when a parseable format is required (e.g. taxonomic_path).
+func (node *TaxNode) FullString(taxonomyCode string) string {
+	if node.HasScientificName() {
+		return fmt.Sprintf("%s:%v [%s]@%s",
+			taxonomyCode,
+			*node.id,
+			node.ScientificName(),
+			node.Rank(),
+		)
+	}
+
+	return fmt.Sprintf("%s:%v",
+		taxonomyCode,
+		*node.id)
+}
+
 // String returns a string representation of the TaxNode, including the taxonomy code,
 // the node ID, and the scientific name. The output format is "taxonomyCode:id [scientificName]".
 //
@@ -42,19 +60,7 @@ func (node *TaxNode) String(taxonomyCode string) string {
 		return *node.id
 	}
 
-	if node.HasScientificName() {
+	return node.FullString(taxonomyCode)
-		return fmt.Sprintf("%s:%v [%s]@%s",
-			taxonomyCode,
-			*node.id,
-			node.ScientificName(),
-			node.Rank(),
-		)
-	}
-
-	return fmt.Sprintf("%s:%v",
-		taxonomyCode,
-		*node.id)
-
 }
 
 // Id returns the unique identifier of the TaxNode.

pkg/obitools/obiconvert/options.go

@@ -24,6 +24,7 @@ var __input_genbank_format__ = false
 var __input_fastq_format__ = false
 var __input_fasta_format__ = false
 var __input_csv_format__ = false
+var __input_json_format__ = false
 
 var __output_in_fasta__ = false
 var __output_in_fastq__ = false
@@ -71,6 +72,9 @@ func InputOptionSet(options *getoptions.GetOpt) {
 	options.BoolVar(&__input_csv_format__, "csv", __input_csv_format__,
 		options.Description("Read data following the CSV format."))
 
+	options.BoolVar(&__input_json_format__, "json", __input_json_format__,
+		options.Description("Read data following the JSON format."))
+
 	options.BoolVar(&__no_ordered_input__, "no-order", __no_ordered_input__,
 		options.Description("When several input files are provided, "+
 			"indicates that there is no order among them."))
@@ -158,6 +162,8 @@ func CLIInputFormat() string {
 		return "genbank"
 	case __input_csv_format__:
 		return "csv"
+	case __input_json_format__:
+		return "json"
 	default:
 		return "guessed"
 	}

pkg/obitools/obiconvert/sequence_reader.go

@@ -73,7 +73,9 @@ func ExpandListOfFiles(check_ext bool, filenames ...string) ([]string, error) {
 						strings.HasSuffix(path, "dat") ||
 						strings.HasSuffix(path, "dat.gz") ||
 						strings.HasSuffix(path, "ecopcr") ||
-						strings.HasSuffix(path, "ecopcr.gz") {
+						strings.HasSuffix(path, "ecopcr.gz") ||
+						strings.HasSuffix(path, "json") ||
+						strings.HasSuffix(path, "json.gz") {
 						log.Debugf("Appending %s file\n", path)
 						list_of_files.Add(path)
 					}
@@ -142,6 +144,8 @@ func CLIReadBioSequences(filenames ...string) (obiiter.IBioSequence, error) {
 			iterator, err = obiformats.ReadFastq(os.Stdin, opts...)
 		case "csv":
 			iterator, err = obiformats.ReadCSV(os.Stdin, opts...)
+		case "json":
+			iterator, err = obiformats.ReadJSON(os.Stdin, opts...)
 		default:
 			iterator, err = obiformats.ReadSequencesFromStdin(opts...)
 		}
@@ -163,6 +167,8 @@ func CLIReadBioSequences(filenames ...string) (obiiter.IBioSequence, error) {
 			reader = obiformats.ReadFastaFromFile
 		case "csv":
 			reader = obiformats.ReadCSVFromFile
+		case "json":
+			reader = obiformats.ReadJSONFromFile
 		case "ecopcr":
 			reader = obiformats.ReadEcoPCRFromFile
 		case "embl":

pkg/obitools/obicsv/obicsv.go

@@ -33,6 +33,7 @@ func CLIWriteSequenceCSV(iterator obiiter.IBioSequence,
 		CSVSequence(CLIPrintSequence()),
 		CSVQuality(CLIPrintQuality()),
 		CSVAutoColumn(CLIAutoColumns()),
+		CSVNAValue(CLINAValue()),
 	)
 
 	csvIter := NewCSVSequenceIterator(iterator, opts...)

pkg/obitools/obicsv/sequence.go

@@ -1,6 +1,7 @@
 package obicsv
 
 import (
+	"fmt"
 	"log"
 	"slices"
 
@@ -67,8 +68,19 @@ func CSVBatchFromSequences(batch obiiter.BioSequenceBatch, opt Options) obiiterc
 
 			if taxon != nil {
 				taxid = taxon.String()
+			} else if ta, ok := sequence.GetAttribute("taxid"); ok {
+				switch tv := ta.(type) {
+				case string:
+					taxid = tv
+				case int:
+					taxid = fmt.Sprintf("%d", tv)
+				case float64:
+					taxid = fmt.Sprintf("%d", int(tv))
+				default:
+					taxid = opt.CSVNAValue()
+				}
 			} else {
-				taxid = sequence.Taxid()
+				taxid = opt.CSVNAValue()
 			}
 
 			record["taxid"] = taxid

pkg/obitools/obidistribute/distribute.go

@@ -46,8 +46,7 @@ func CLIDistributeSequence(sequences obiiter.IBioSequence) {
 		formater = obiformats.WriteSequencesToFile
 	}
 
-	dispatcher := sequences.Distribute(CLISequenceClassifier(),
+	dispatcher := sequences.Distribute(CLISequenceClassifier())
-		obidefault.BatchSize())
 
 	obiformats.WriterDispatcher(CLIFileNamePattern(),
 		dispatcher, formater, opts...,

pkg/obitools/obilandmark/obilandmark.go

@@ -170,7 +170,7 @@ func CLISelectLandmarkSequences(iterator obiiter.IBioSequence) obiiter.IBioSeque
 		for i, seq := range library {
 			taxon := seq.Taxon(taxo)
 			if taxon == nil {
-				log.Fatal("%s: Cannot identify taxid %s in %s", seq.Id(), seq.Taxid(), taxo.Name())
+				log.Fatalf("%s: Cannot identify taxid %s in %s", seq.Id(), seq.Taxid(), taxo.Name())
 			}
 			taxa.Set(i, taxon)
 		}

pkg/obitools/obimultiplex/demultiplex.go

@@ -15,7 +15,6 @@ func IExtractBarcode(iterator obiiter.IBioSequence) (obiiter.IBioSequence, error
 	opts := make([]obingslibrary.WithOption, 0, 10)
 
 	opts = append(opts,
-		obingslibrary.OptionAllowedMismatches(CLIAllowedMismatch()),
 		obingslibrary.OptionAllowedIndel(CLIAllowsIndel()),
 		obingslibrary.OptionUnidentified(CLIUnidentifiedFileName()),
 		obingslibrary.OptionDiscardErrors(!CLIConservedErrors()),
@@ -23,6 +22,14 @@ func IExtractBarcode(iterator obiiter.IBioSequence) (obiiter.IBioSequence, error
 		obingslibrary.OptionBatchSize(obidefault.BatchSize()),
 	)
 
+	// Only propagate the CLI --allowed-mismatches value if the user
+	// explicitly set it: otherwise the per-primer values defined in
+	// the NGSFilter config file (@primer_mismatches, @forward_mismatches,
+	// @reverse_mismatches) must be preserved.
+	if CLIAllowedMismatchIsSet() {
+		opts = append(opts, obingslibrary.OptionAllowedMismatches(CLIAllowedMismatch()))
+	}
+
 	ngsfilter, err := CLINGSFIlter()
 	if err != nil {
 		log.Fatalf("%v", err)

pkg/obitools/obimultiplex/options.go

@@ -18,6 +18,7 @@ var _UnidentifiedFile = ""
 var _AllowedMismatch = 2
 var _AllowsIndel = false
 var _ConservedError = false
+var _optionsParser *getoptions.GetOpt
 
 // PCROptionSet defines every options related to a simulated PCR.
 //
@@ -29,6 +30,8 @@ var _ConservedError = false
 // - option : is a pointer to a getoptions.GetOpt instance normaly
 // produced by the
 func MultiplexOptionSet(options *getoptions.GetOpt) {
+	_optionsParser = options
+
 	options.StringVar(&_NGSFilterFile, "tag-list", _NGSFilterFile,
 		options.Alias("s"),
 		options.Description("File name of the NGSFilter file describing PCRs."))
@@ -62,6 +65,15 @@ func CLIAllowedMismatch() int {
 	return _AllowedMismatch
 }
 
+// CLIAllowedMismatchIsSet returns true if the user explicitly
+// specified --allowed-mismatches on the command line, as opposed
+// to relying on its default value. This allows per-primer mismatch
+// settings from the NGSFilter config file to take precedence unless
+// the user explicitly overrides them from the CLI.
+func CLIAllowedMismatchIsSet() bool {
+	return _optionsParser != nil && _optionsParser.Called("allowed-mismatches")
+}
+
 func CLIAllowsIndel() bool {
 	return _AllowsIndel
 }

pkg/obitools/obipairing/options.go

@@ -21,12 +21,10 @@ func PairingOptionSet(options *getoptions.GetOpt) {
 	options.StringVar(&_ForwardFile, "forward-reads", "",
 		options.Alias("F"),
 		options.ArgName("FILENAME_F"),
-		options.Required("You must provide at a forward file"),
 		options.Description("The file names containing the forward reads"))
 	options.StringVar(&_ReverseFile, "reverse-reads", "",
 		options.Alias("R"),
 		options.ArgName("FILENAME_R"),
-		options.Required("You must provide a reverse file"),
 		options.Description("The file names containing the reverse reads"))
 	options.IntVar(&_Delta, "delta", _Delta,
 		options.Alias("D"),
@@ -72,6 +70,10 @@ func CLIPairedSequence() (obiiter.IBioSequence, error) {
 	return paired, nil
 }
 
+func CLIHasPairedFiles() bool {
+	return _ForwardFile != "" && _ReverseFile != ""
+}
+
 func CLIDelta() int {
 	return _Delta
 }

pkg/obitools/obisummary/obisummary.go

@@ -99,6 +99,17 @@ func (data1 *DataSummary) Add(data2 *DataSummary) *DataSummary {
 	rep.sample_singletons = sumUpdateIntMap(data1.sample_singletons, data2.sample_singletons)
 	rep.sample_obiclean_bad = sumUpdateIntMap(data1.sample_obiclean_bad, data2.sample_obiclean_bad)
 
+	for k, m1 := range data1.map_summaries {
+		rep.map_summaries[k] = m1
+	}
+	for k, m2 := range data2.map_summaries {
+		if m1, ok := rep.map_summaries[k]; ok {
+			rep.map_summaries[k] = sumUpdateIntMap(m1, m2)
+		} else {
+			rep.map_summaries[k] = m2
+		}
+	}
+
 	return rep
 }
 
@@ -163,8 +174,9 @@ func ISummary(iterator obiiter.IBioSequence, summarise []string) map[string]inte
 	summaries := make([]*DataSummary, nproc)
 
 	for n := 0; n < nproc; n++ {
+		summaries[n] = NewDataSummary()
 		for _, v := range summarise {
-			summaries[n].map_summaries[v] = make(map[string]int, 0)
+			summaries[n].map_summaries[v] = make(map[string]int)
 		}
 	}
 
@@ -174,6 +186,11 @@ func ISummary(iterator obiiter.IBioSequence, summarise []string) map[string]inte
 			batch := iseq.Get()
 			for _, seq := range batch.Slice() {
 				summary.Update(seq)
+				for _, attr := range summarise {
+					if m, ok := seq.GetIntMap(attr); ok {
+						summary.map_summaries[attr] = sumUpdateIntMap(summary.map_summaries[attr], m)
+					}
+				}
 			}
 		}
 		waiter.Done()
@@ -181,11 +198,9 @@ func ISummary(iterator obiiter.IBioSequence, summarise []string) map[string]inte
 
 	waiter.Add(nproc)
 
-	summaries[0] = NewDataSummary()
 	go ff(iterator, summaries[0])
 
 	for i := 1; i < nproc; i++ {
-		summaries[i] = NewDataSummary()
 		go ff(iterator.Split(), summaries[i])
 	}
 
@@ -246,5 +261,14 @@ func ISummary(iterator obiiter.IBioSequence, summarise []string) map[string]inte
 			}
 		}
 	}
+
+	if len(rep.map_summaries) > 0 {
+		mapDict := make(map[string]interface{}, len(rep.map_summaries))
+		for attr, counts := range rep.map_summaries {
+			mapDict[attr] = counts
+		}
+		dict["map_summaries"] = mapDict
+	}
+
 	return dict
 }

pkg/obitools/obitagpcr/pcrtag.go

@@ -55,7 +55,7 @@ func IPCRTagPESequencesBatch(iterator obiiter.IBioSequence,
 		ngsfilter.SetAllowsIndels(true)
 	}
 
-	if obimultiplex.CLIAllowedMismatch() > 0 {
+	if obimultiplex.CLIAllowedMismatchIsSet() {
 		ngsfilter.SetAllowedMismatches(obimultiplex.CLIAllowedMismatch())
 	}
 
@@ -114,10 +114,10 @@ func IPCRTagPESequencesBatch(iterator obiiter.IBioSequence,
 					aanot["obimultiplex_direction"] = direction
 
 					aanot["obimultiplex_forward_match"] = forward_match
-					aanot["obimultiplex_forward_mismatches"] = forward_mismatches
+					aanot["obimultiplex_forward_error"] = forward_mismatches
 
 					aanot["obimultiplex_reverse_match"] = reverse_match
-					aanot["obimultiplex_reverse_mismatches"] = reverse_mismatches
+					aanot["obimultiplex_reverse_error"] = reverse_mismatches
 
 					aanot["sample"] = sample
 					aanot["experiment"] = experiment
@@ -125,10 +125,10 @@ func IPCRTagPESequencesBatch(iterator obiiter.IBioSequence,
 					banot["obimultiplex_direction"] = direction
 
 					banot["obimultiplex_forward_match"] = forward_match
-					banot["obimultiplex_forward_mismatches"] = forward_mismatches
+					banot["obimultiplex_forward_error"] = forward_mismatches
 
 					banot["obimultiplex_reverse_match"] = reverse_match
-					banot["obimultiplex_reverse_mismatches"] = reverse_mismatches
+					banot["obimultiplex_reverse_error"] = reverse_mismatches
 
 					banot["sample"] = sample
 					banot["experiment"] = experiment

pkg/obiutils/cast_interface.go

@@ -276,6 +276,44 @@ func InterfaceToStringMap(i interface{}) (val map[string]string, err error) {
 	return
 }
 
+func InterfaceToMapOfIntSlice(i interface{}) (val map[string][]int, err error) {
+	err = nil
+	switch m := i.(type) {
+	case map[string][]int:
+		val = m
+	case map[string]interface{}:
+		val = make(map[string][]int, len(m))
+		for k, v := range m {
+			val[k], err = InterfaceToIntSlice(v)
+			if err != nil {
+				return
+			}
+		}
+	default:
+		err = &NotAMapInt{"value attribute cannot be casted to a map[string][]int"}
+	}
+	return
+}
+
+func InterfaceToMapOfStringSlice(i interface{}) (val map[string][]string, err error) {
+	err = nil
+	switch m := i.(type) {
+	case map[string][]string:
+		val = m
+	case map[string]interface{}:
+		val = make(map[string][]string, len(m))
+		for k, v := range m {
+			val[k], err = InterfaceToStringSlice(v)
+			if err != nil {
+				return
+			}
+		}
+	default:
+		err = &NotAMapInt{"value attribute cannot be casted to a map[string][]string"}
+	}
+	return
+}
+
 func InterfaceToStringSlice(i interface{}) (val []string, err error) {
 	err = nil
 

pkg/obiutils/mimetypes.go

@@ -102,6 +102,11 @@ func RegisterOBIMimeType() {
 			return ok
 		}
 
+		jsonDetector := func(raw []byte, limit uint32) bool {
+			raw = bytes.TrimLeft(raw, " \t\r\n")
+			return len(raw) > 0 && (raw[0] == '[' || raw[0] == '{')
+		}
+
 		mimetype.Lookup("text/plain").Extend(fastaDetector, "text/fasta", ".fasta")
 		mimetype.Lookup("text/plain").Extend(fastqDetector, "text/fastq", ".fastq")
 		mimetype.Lookup("text/plain").Extend(ecoPCR2Detector, "text/ecopcr2", ".ecopcr")
@@ -115,6 +120,7 @@ func RegisterOBIMimeType() {
 		mimetype.Lookup("application/octet-stream").Extend(genbankDetector, "text/genbank", ".seq")
 		mimetype.Lookup("application/octet-stream").Extend(emblDetector, "text/embl", ".dat")
 		mimetype.Lookup("application/octet-stream").Extend(csv, "text/csv", ".csv")
+		mimetype.Lookup("application/octet-stream").Extend(jsonDetector, "application/json", ".json")
 	}
 	__obimimetype_registred__ = true
 }

pkg/obiutils/minmax.go

@@ -34,6 +34,26 @@ func MinMaxSlice[T constraints.Ordered](vec []T) (min, max T) {
 	return
 }
 
+func FilterMinSlice[T constraints.Ordered](vec []T, minimum T) []T {
+	result := make([]T, 0, len(vec))
+	for _, v := range vec {
+		if v >= minimum {
+			result = append(result, v)
+		}
+	}
+	return result
+}
+
+func FilterMaxSlice[T constraints.Ordered](vec []T, maximum T) []T {
+	result := make([]T, 0, len(vec))
+	for _, v := range vec {
+		if v <= maximum {
+			result = append(result, v)
+		}
+	}
+	return result
+}
+
 func MaxMap[K comparable, T constraints.Ordered](values map[K]T) (K, T, error) {
 	var maxKey K
 	var maxValue T
@@ -73,6 +93,46 @@ func MinMap[K comparable, T constraints.Ordered](values map[K]T) (K, T, error) {
 	return minKey, minValue, nil
 }
 
+func FilterMinMap[K comparable, T constraints.Ordered](values map[K]T, minimum T) map[K]T {
+	result := make(map[K]T)
+	for k, v := range values {
+		if v >= minimum {
+			result[k] = v
+		}
+	}
+	return result
+}
+
+func FilterMaxMap[K comparable, T constraints.Ordered](values map[K]T, maximum T) map[K]T {
+	result := make(map[K]T)
+	for k, v := range values {
+		if v <= maximum {
+			result[k] = v
+		}
+	}
+	return result
+}
+
+func SaturatingSubSlice[T Numeric](vec []T, sub T) []T {
+	result := make([]T, len(vec))
+	for i, v := range vec {
+		if v > sub {
+			result[i] = v - sub
+		}
+	}
+	return result
+}
+
+func SaturatingSubMap[K comparable, T Numeric](values map[K]T, sub T) map[K]T {
+	result := make(map[K]T)
+	for k, v := range values {
+		if v > sub {
+			result[k] = v - sub
+		}
+	}
+	return result
+}
+
 // Min returns the smallest element in a slice/array or map,
 // or the value itself if data is a single comparable value.
 // Returns an error if the container is empty or the type is unsupported.
@@ -135,11 +195,121 @@ func Max(data interface{}) (interface{}, error) {
 	}
 }
 
+func FilterMin(data interface{}, minimum interface{}) (interface{}, error) {
+	v := reflect.ValueOf(data)
+	switch v.Kind() {
+	case reflect.Slice, reflect.Array:
+		if v.Len() == 0 {
+			return nil, errors.New("empty slice or array")
+		}
+		return filterMinFromIterable(v, minimum)
+	case reflect.Map:
+		if v.Len() == 0 {
+			return nil, errors.New("empty map")
+		}
+		return filterMinFromMap(v, minimum)
+	default:
+		if !isOrderedKind(v.Kind()) {
+			return nil, fmt.Errorf("unsupported type: %s", v.Kind())
+		}
+		return data, nil
+	}
+}
+
+func FilterMax(data interface{}, maximum interface{}) (interface{}, error) {
+	v := reflect.ValueOf(data)
+	switch v.Kind() {
+	case reflect.Slice, reflect.Array:
+		if v.Len() == 0 {
+			return nil, errors.New("empty slice or array")
+		}
+		return filterMaxFromIterable(v, maximum)
+	case reflect.Map:
+		if v.Len() == 0 {
+			return nil, errors.New("empty map")
+		}
+		return filterMaxFromMap(v, maximum)
+	default:
+		if !isOrderedKind(v.Kind()) {
+			return nil, fmt.Errorf("unsupported type: %s", v.Kind())
+		}
+		return data, nil
+	}
+}
+
+func SaturatingSub(data interface{}, sub interface{}) (interface{}, error) {
+	v := reflect.ValueOf(data)
+	switch v.Kind() {
+	case reflect.Slice, reflect.Array:
+		return saturatingSubFromIterable(v, sub)
+	case reflect.Map:
+		return saturatingSubFromMap(v, sub)
+	default:
+		if !isNumericKind(v.Kind()) {
+			return nil, fmt.Errorf("unsupported type: %s", v.Kind())
+		}
+		r, err := saturatingSubValues(v, reflect.ValueOf(sub))
+		if err != nil {
+			return nil, err
+		}
+		return r.Interface(), nil
+	}
+}
+
+func saturatingSubFromIterable(v reflect.Value, sub interface{}) (interface{}, error) {
+	subVal := reflect.ValueOf(sub)
+	result := reflect.MakeSlice(v.Type(), v.Len(), v.Len())
+	for i := 0; i < v.Len(); i++ {
+		r, err := saturatingSubValues(v.Index(i), subVal)
+		if err != nil {
+			return nil, err
+		}
+		result.Index(i).Set(r)
+	}
+	return result.Interface(), nil
+}
+
+func saturatingSubFromMap(v reflect.Value, sub interface{}) (interface{}, error) {
+	subVal := reflect.ValueOf(sub)
+	result := reflect.MakeMap(v.Type())
+	for _, key := range v.MapKeys() {
+		r, err := saturatingSubValues(v.MapIndex(key), subVal)
+		if err != nil {
+			return nil, err
+		}
+		if !r.IsZero() {
+			result.SetMapIndex(key, r)
+		}
+	}
+	return result.Interface(), nil
+}
+
+func saturatingSubValues(a, b reflect.Value) (reflect.Value, error) {
+	result := reflect.New(a.Type()).Elem()
+	switch a.Kind() {
+	case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64:
+		if av, bv := a.Int(), b.Int(); av > bv {
+			result.SetInt(av - bv)
+		}
+	case reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64:
+		if av, bv := a.Uint(), b.Uint(); av > bv {
+			result.SetUint(av - bv)
+		}
+	case reflect.Float32, reflect.Float64:
+		if av, bv := a.Float(), b.Float(); av > bv {
+			result.SetFloat(av - bv)
+		}
+	default:
+		return reflect.Value{}, fmt.Errorf("unsupported type for saturating subtraction: %s", a.Kind())
+	}
+	return result, nil
+}
+
 // maxFromIterable scans a slice/array to find the maximum.
 func maxFromIterable(v reflect.Value) (interface{}, error) {
 	var best reflect.Value
 	for i := 0; i < v.Len(); i++ {
-		elem := v.Index(i)
+		elem := unwrapInterface(v.Index(i))
 		if !isOrderedKind(elem.Kind()) {
 			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
 		}
@@ -154,7 +324,7 @@ func maxFromIterable(v reflect.Value) (interface{}, error) {
 func minFromIterable(v reflect.Value) (interface{}, error) {
 	var minVal reflect.Value
 	for i := 0; i < v.Len(); i++ {
-		elem := v.Index(i)
+		elem := unwrapInterface(v.Index(i))
 		if !isOrderedKind(elem.Kind()) {
 			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
 		}
@@ -165,12 +335,182 @@ func minFromIterable(v reflect.Value) (interface{}, error) {
 	return minVal.Interface(), nil
 }
 
+func filterMinFromIterable(v reflect.Value, minimum interface{}) (interface{}, error) {
+	minVal := reflect.ValueOf(minimum)
+	result := reflect.MakeSlice(v.Type(), 0, v.Len())
+	for i := 0; i < v.Len(); i++ {
+		elem := unwrapInterface(v.Index(i))
+		if !isOrderedKind(elem.Kind()) {
+			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
+		}
+		if !less(elem, minVal) { // elem >= minimum
+			result = reflect.Append(result, elem)
+		}
+	}
+	return result.Interface(), nil
+}
+
+func filterMaxFromIterable(v reflect.Value, maximum interface{}) (interface{}, error) {
+	maxVal := reflect.ValueOf(maximum)
+	result := reflect.MakeSlice(v.Type(), 0, v.Len())
+	for i := 0; i < v.Len(); i++ {
+		elem := unwrapInterface(v.Index(i))
+		if !isOrderedKind(elem.Kind()) {
+			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
+		}
+		if !greater(elem, maxVal) { // elem <= maximum
+			result = reflect.Append(result, elem)
+		}
+	}
+	return result.Interface(), nil
+}
+
+// whichMaxFromIterable returns the index of the maximum element in a slice/array.
+func whichMaxFromIterable(v reflect.Value) (int, error) {
+	var best reflect.Value
+	bestIdx := 0
+	for i := 0; i < v.Len(); i++ {
+		elem := unwrapInterface(v.Index(i))
+		if !isOrderedKind(elem.Kind()) {
+			return 0, fmt.Errorf("unsupported element type: %s", elem.Kind())
+		}
+		if i == 0 || greater(elem, best) {
+			best = elem
+			bestIdx = i
+		}
+	}
+	return bestIdx, nil
+}
+
+// whichMinFromIterable returns the index of the minimum element in a slice/array.
+func whichMinFromIterable(v reflect.Value) (int, error) {
+	var minVal reflect.Value
+	minIdx := 0
+	for i := 0; i < v.Len(); i++ {
+		elem := unwrapInterface(v.Index(i))
+		if !isOrderedKind(elem.Kind()) {
+			return 0, fmt.Errorf("unsupported element type: %s", elem.Kind())
+		}
+		if i == 0 || less(elem, minVal) {
+			minVal = elem
+			minIdx = i
+		}
+	}
+	return minIdx, nil
+}
+
+// whichMaxFromMap returns the key associated with the maximum value in a map.
+func whichMaxFromMap(v reflect.Value) (interface{}, error) {
+	var best reflect.Value
+	var bestKey reflect.Value
+	first := true
+	for _, key := range v.MapKeys() {
+		elem := unwrapInterface(v.MapIndex(key))
+		if !isOrderedKind(elem.Kind()) {
+			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
+		}
+		if first || greater(elem, best) {
+			best = elem
+			bestKey = key
+			first = false
+		}
+	}
+	return bestKey.Interface(), nil
+}
+
+// whichMinFromMap returns the key associated with the minimum value in a map.
+func whichMinFromMap(v reflect.Value) (interface{}, error) {
+	var minVal reflect.Value
+	var minKey reflect.Value
+	first := true
+	for _, key := range v.MapKeys() {
+		elem := unwrapInterface(v.MapIndex(key))
+		if !isOrderedKind(elem.Kind()) {
+			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
+		}
+		if first || less(elem, minVal) {
+			minVal = elem
+			minKey = key
+			first = false
+		}
+	}
+	return minKey.Interface(), nil
+}
+
+// WhichMax returns the key (for a map) or index (for a slice/array) of the maximum value.
+func WhichMax(data interface{}) (interface{}, error) {
+	v := reflect.ValueOf(data)
+	switch v.Kind() {
+	case reflect.Slice, reflect.Array:
+		if v.Len() == 0 {
+			return nil, errors.New("empty slice or array")
+		}
+		return whichMaxFromIterable(v)
+	case reflect.Map:
+		if v.Len() == 0 {
+			return nil, errors.New("empty map")
+		}
+		return whichMaxFromMap(v)
+	default:
+		return nil, fmt.Errorf("unsupported type: %s", v.Kind())
+	}
+}
+
+// WhichMin returns the key (for a map) or index (for a slice/array) of the minimum value.
+func WhichMin(data interface{}) (interface{}, error) {
+	v := reflect.ValueOf(data)
+	switch v.Kind() {
+	case reflect.Slice, reflect.Array:
+		if v.Len() == 0 {
+			return nil, errors.New("empty slice or array")
+		}
+		return whichMinFromIterable(v)
+	case reflect.Map:
+		if v.Len() == 0 {
+			return nil, errors.New("empty map")
+		}
+		return whichMinFromMap(v)
+	default:
+		return nil, fmt.Errorf("unsupported type: %s", v.Kind())
+	}
+}
+
+func filterMinFromMap(v reflect.Value, minimum interface{}) (interface{}, error) {
+	minVal := reflect.ValueOf(minimum)
+	result := reflect.MakeMap(v.Type())
+	for _, key := range v.MapKeys() {
+		elem := unwrapInterface(v.MapIndex(key))
+		if !isOrderedKind(elem.Kind()) {
+			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
+		}
+		if !less(elem, minVal) { // elem >= minimum
+			result.SetMapIndex(key, elem)
+		}
+	}
+	return result.Interface(), nil
+}
+
+func filterMaxFromMap(v reflect.Value, maximum interface{}) (interface{}, error) {
+	maxVal := reflect.ValueOf(maximum)
+	result := reflect.MakeMap(v.Type())
+	for _, key := range v.MapKeys() {
+		elem := unwrapInterface(v.MapIndex(key))
+		if !isOrderedKind(elem.Kind()) {
+			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
+		}
+		if !greater(elem, maxVal) { // elem <= maximum
+			result.SetMapIndex(key, elem)
+		}
+	}
+	return result.Interface(), nil
+}
+
 // maxFromMap scans map values to find the maximum.
 func maxFromMap(v reflect.Value) (interface{}, error) {
 	var best reflect.Value
 	first := true
 	for _, key := range v.MapKeys() {
-		elem := v.MapIndex(key)
+		elem := unwrapInterface(v.MapIndex(key))
 		if !isOrderedKind(elem.Kind()) {
 			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
 		}
@@ -187,7 +527,7 @@ func minFromMap(v reflect.Value) (interface{}, error) {
 	var minVal reflect.Value
 	first := true
 	for _, key := range v.MapKeys() {
-		elem := v.MapIndex(key)
+		elem := unwrapInterface(v.MapIndex(key))
 		if !isOrderedKind(elem.Kind()) {
 			return nil, fmt.Errorf("unsupported element type: %s", elem.Kind())
 		}
@@ -199,6 +539,27 @@ func minFromMap(v reflect.Value) (interface{}, error) {
 	return minVal.Interface(), nil
 }
 
+func isNumericKind(k reflect.Kind) bool {
+	switch k {
+	case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64,
+		reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64,
+		reflect.Float32, reflect.Float64:
+		return true
+	default:
+		return false
+	}
+}
+
+// unwrapInterface returns v.Elem() when v holds an interface value, otherwise v unchanged.
+// This is necessary when iterating map[string]interface{} or []interface{} via reflection:
+// the element Kind is reflect.Interface, not the underlying concrete type.
+func unwrapInterface(v reflect.Value) reflect.Value {
+	if v.Kind() == reflect.Interface {
+		return v.Elem()
+	}
+	return v
+}
+
 // isOrderedKind reports whether k supports comparison ordering.
 func isOrderedKind(k reflect.Kind) bool {
 	switch k {

prompt_documentation_globale.md

@@ -0,0 +1,302 @@
+# Objective
+
+Fully document OBITools (version 4, written in Go) in English, using a 4‑phase incremental pipeline.
+
+You **MUST** use the available MCP servers:
+
+- `cclsp` – exact definitions, references, diagnostics  
+- `jcodemunch` – code indexing, symbol extraction  
+- `treesitter` – AST and CLI parsing  
+- `context7` – external documentation  
+
+All tool calls must follow the exact API described in the MCP server documentation. If a required tool is unavailable, you **MUST** log the error and stop execution.
+
+### Tool call format (CRITICAL)
+
+Tool calls **MUST** use this exact XML format — no spaces inside the angle brackets:
+
+```
+<function=tool_name>
+{"param": "value"}
+</function>
+```
+
+**FORBIDDEN** — these variants will cause parse errors and must NEVER be used:
+- `< function=tool_name >` (spaces around the tag name)
+- `< function = tool_name>` (spaces around `=`)
+- `<function = tool_name>` (space before `=`)
+
+The opening tag is `<function=tool_name>` with **zero spaces** inside `<` and `>`.
+
+---
+
+# Global rules
+
+** You are not allowed to read twice the same file in a row. **
+
+## Language
+
+- All generated documentation **MUST** be in English.  
+- If an existing documentation file is in French:  
+  1. Translate it to English  
+  2. Save the original as `.fr.md` **before** overwriting  
+  3. Write the new English version  
+
+---
+
+## Execution mode (STRICT)
+
+You are operating in **STRICT TOOL MODE**:
+
+- If a file must be written, you **MUST** use the `Shell` tool.
+- You **MUST NOT** read entire directory listings into memory.
+- You **MUST** work with **one item at a time** using a simple text file as a task queue.
+
+### Reading files before writing
+
+- **Before writing to an existing documentation file**, you must first read it using the `Read` tool.
+- **When documenting a single Go source file**, you only need to read that one file (plus up to 4-5 helper files if needed for context).
+- Do NOT read the entire codebase - only what is necessary to document the current file.
+
+---
+
+### Rules
+
+- Always write the **full** file (no partial updates).  
+- Paths are relative to the project root; directories are created implicitly.  
+- Content must be valid UTF‑8; use `\n` line endings.  
+- Do **not** wrap content in backticks.  
+
+---
+
+## Progress tracking: task queue files
+
+We use **line‑oriented task files** to avoid loading large lists into memory. Each phase has its own task file:
+
+- `docs/todo/phase1.txt` – list of Go files (one per line) to document.  
+- `docs/todo/phase1bis.txt` – same list, but after phase1 is done.  
+- `docs/todo/phase2.txt` – list of packages.  
+- `docs/todo/phase3.txt` – list of tools.
+
+**How it works:**
+
+1. At the start of a phase, if the task file does not exist, it is created by scanning the codebase once (Phase 0 or Phase X init).  
+2. **Each run of the LLM processes only the first line of the task file.**  
+3. After processing the item (success or permanent failure), the line is removed from the task file.  
+   - On success, the line is deleted (no extra sentinel file needed).  
+   - On transient failure (retry < 3), we keep the line but increment a retry counter stored in a separate file.  
+   - On permanent failure (retry ≥ 3), we move the line to a `failed.txt` file and log the error.  
+4. The LLM then exits (or continues if the task file is still non‑empty, but it must never load more than one line).
+
+This way, the LLM’s context never holds more than a single task at a time.
+
+### Retry mechanism
+
+For each item (e.g., `internal/align/align.go`), we maintain a retry counter in:
+
+- `docs/retry/phase1/internal/align/align.go.count`
+
+If the file does not exist, retries = 0.  
+Each time processing fails, we increment the counter (write the new number).  
+If after increment the counter < 3, we keep the line in the task file.  
+If counter reaches 3, we **remove the line from the task file**, add it to `docs/failed/phase1/internal/align/align.go.failed` (just a marker), and log the error.
+
+---
+
+## Documentation quality requirements (CRITICAL)
+
+Documentation MUST NOT be superficial. For each documented element (file, function, struct, package):
+
+###  You MUST explain:
+
+- what it does
+- why it exists (context, problem solved)
+- how it is used
+- assumptions and preconditions
+- possible edge cases
+
+### Forbidden patterns
+
+- Vague phrases like “This function handles…”, “Utility for…”, “Helper function…”.
+- Generic descriptions that could apply to any project.
+
+### Required content per element type
+
+- Functions:
+  - Purpose
+  - Parameter meaning
+  - Return values
+  - Notable behaviour (panic conditions, side effects, concurrency)
+- Structs:
+  - Role in the system
+  - Meaning of key fields
+- Files:
+  - Role within the package
+  - Interactions with other files
+
+### Anti‑generic rule
+
+If the description could apply to any project, it is INVALID. You MUST include domain‑specific context (bioinformatics, sequence processing, etc.) and concrete behaviour.
+
+### Quality validation
+
+Before marking an item as done (i.e., creating the .done sentinel), you MUST perform a self‑validation:
+
+- Check that all required sections are present.
+- Verify that no forbidden patterns remain.
+
+If validation fails, increment the retry counter and keep the item pending.
+
+
+---
+
+# Directory structure
+
+```
+docs/
+  todo/                          # task queues
+    phase1.txt
+    phase1bis.txt
+    phase2.txt
+    phase3.txt
+  retry/                         # retry counters
+    phase1/                      # mirrors file structure
+      internal/align/align.go.count
+    phase1bis/
+    phase2/
+    phase3/
+  failed/                        # permanent failure markers
+    phase1/
+      internal/align/align.go.failed
+    phase1bis/
+    phase2/
+    phase3/
+  phase1/                        # actual documentation
+    <relative_path>/<file>.go.md
+  phase2/
+    <package>.md
+  phase3/
+    <tool>.md
+  error.log
+```
+
+---
+
+# Phase 0: Initialization
+
+1. Ensure required directories exist: `docs/todo`, `docs/retry`, `docs/failed`, `docs/phase1`, `docs/phase2`, `docs/phase3`.  
+2. **If `docs/todo/phase1.txt` does not exist**:  
+   - Use `find pkg -name "*.go" ! -name "*_test.go" ! -path "*/cmd/*"` to list all Go files (excluding tests and main.go).  
+   - Write the list (one relative path per line, e.g., `internal/align/align.go`) to `docs/todo/phase1.txt`.  
+3. Do the same for phase2 and phase3 later when those phases start.  
+4. **No other state is stored.**
+
+---
+
+# Phase 1: File documentation
+
+**Processing rule:**
+- Read the **first line** of `docs/todo/phase1.txt` (using `head -n 1`).  
+- If the file is empty, Phase 1 is complete → proceed to Phase 1bis initialization.  
+- Otherwise, process that single file.
+
+**Processing a file:**
+
+1. Let `relpath` be the line content (e.g., `internal/align/align.go`).  
+2. Check if a permanent failure marker exists at `docs/failed/phase1/${relpath}.failed`. If yes, remove the line from the task file and skip (line will be deleted). 
+3. If the documentation file `docs/phase1/${relpath}.go.md` exists go directly to its validation (step 6).
+4. Otherwise, generate documentation for that file (using MCP tools as before).  
+5. Write the documentation to `docs/phase1/${relpath}.go.md`.  
+6. Validate quality.  
+7. If validation succeeds:  
+   - Remove the line from the task file.  
+   - Remove any retry counter file for this item.  
+   - (No sentinel needed; the removal from todo indicates completion.)  
+8. If validation fails:  
+   - Increment retry counter:  
+     - If `docs/retry/phase1/${relpath}.count` does not exist, set to 1.  
+     - Else read it, add 1, write back.  
+   - If new counter >= 3:  
+     - Remove line from task file.  
+     - Create `docs/failed/phase1/${relpath}.failed`.  
+     - Log error.  
+   - If new counter < 3:  
+     - Keep the line in the task file (do nothing, it stays as first line for next run).  
+9. **Exit** (or stop if this was a single run). The next invocation will read the first line again (same if retry, or next if removed).
+
+**Important:**  
+- Do **not** read more than one line.  
+- Do **not** attempt to process multiple items in one run.  
+- The LLM should finish after handling one item.
+
+---
+
+# Phase 1bis: Review and harmonization
+
+When Phase 1 is complete (i.e., `docs/todo/phase1.txt` empty), we initialize `docs/todo/phase1bis.txt` with the same list of files (the ones that succeeded).  
+But note: we need to know which files were successfully documented. Since we removed lines from `phase1.txt` on success, we need a record. The simplest is to reuse the same list but we can generate it by listing the existing `.go.md` files in `docs/phase1/` (since every successful file has a `.go.md`).  
+Thus, Phase 1bis initialization:
+
+- If `docs/todo/phase1bis.txt` does not exist, create it by listing all `.go.md` files under `docs/phase1/`, stripping the `docs/phase1/` prefix and the `.go.md` suffix, and writing the relative path (same format as phase1).  
+
+Then processing is identical to Phase 1, but using `docs/todo/phase1bis.txt` and output is overwriting the same `.go.md` files (with improvements). Retry counters go in `docs/retry/phase1bis/`.
+
+---
+
+# Phase 2: Package documentation
+
+When Phase 1bis is complete (`docs/todo/phase1bis.txt` empty), initialize `docs/todo/phase2.txt`:
+
+- List all packages: unique directories under `pkg/` that contain at least one `.go` file and are not tools.  
+- Write each package identifier (e.g., `align`, `internal/align`) as a line.
+
+Processing: read first line, generate `docs/phase2/<package>.md`, validate, remove line on success, retry logic in `docs/retry/phase2/`.
+
+---
+
+# Phase 3: Tool documentation
+
+When Phase 2 complete, initialize `docs/todo/phase3.txt`:
+
+- List all directories under `cmd/` that contain a `main.go`. Write each tool name as a line.
+
+Processing: read first line, generate `docs/phase3/<tool>.md`, validate, remove line on success, retry logic in `docs/retry/phase3/`.
+
+---
+
+# Finalization
+
+When all task files are empty and no pending phases, generate `docs/README.md` by:
+
+- Listing all package docs (files in `docs/phase2/`) and linking.  
+- Listing all tool docs (files in `docs/phase3/`) and linking.  
+
+Write using `Shell`.
+
+---
+
+# Execution flow summary
+
+1. **Phase 0**: Create directories and initial `todo/phase1.txt` if missing. Exit.  
+2. **Phase 1**:  
+   - If `todo/phase1.txt` exists and non‑empty → process first line.  
+   - Else → move to Phase 1bis initialization.  
+3. **Phase 1bis**:  
+   - If `todo/phase1bis.txt` does not exist → create from successful phase1 docs.  
+   - If non‑empty → process first line.  
+   - Else → move to Phase 2 initialization.  
+4. **Phase 2**: similar.  
+5. **Phase 3**: similar.  
+6. **Finalization**: generate README.
+
+The LLM should be invoked repeatedly (e.g., by a scheduler) until all phases are done. Each invocation processes exactly one item.
+
+---
+
+# Important reminders
+
+- Always call `Shell` to write files; never output content in plain text.  
+- Validate quality before removing a line from the task file.  
+- Log all failures to `docs/error.log` in JSON lines format.  
+- If any MCP tool fails, treat as failure and increment retry counter.  
+- Never read more than one line from a task file in a single run.

scripts/find_setattribute.go

@@ -0,0 +1,222 @@
+//go:build ignore
+
+package main
+
+import (
+	"encoding/json"
+	"fmt"
+	"os"
+	"os/exec"
+	"path/filepath"
+	"regexp"
+	"sort"
+	"strconv"
+	"strings"
+)
+
+type Reference struct {
+	File     string `json:"file"`
+	Line     int    `json:"line"`
+	Column   int    `json:"column"`
+	Key      string `json:"key"`
+	Function string `json:"function"`
+	Context  string `json:"context"`
+}
+
+type Result struct {
+	Method     string      `json:"method"`
+	Signature  string      `json:"signature"`
+	Definition string      `json:"definition"`
+	References []Reference `json:"references"`
+	Total      int         `json:"total"`
+}
+
+var basePath = "/Users/coissac/Sync/travail/__MOI__/GO/obitools4"
+
+func main() {
+	cmd := exec.Command("rg", "-n", `\.SetAttribute\(`, basePath+"/pkg", "--type", "go")
+	output, err := cmd.Output()
+	if err != nil {
+		fmt.Fprintf(os.Stderr, "Error running rg: %v\n", err)
+		os.Exit(1)
+	}
+
+	lines := strings.Split(string(output), "\n")
+	lineRe := regexp.MustCompile(`^(.+?):(\d+):\s*(.+)$`)
+	keyRe := regexp.MustCompile(`SetAttribute\("([^"]+)"`)
+	templateKeyRe := regexp.MustCompile(`SetAttribute\("([^"]+)[^"]*"\s*,`)
+
+	var refs []Reference
+	seen := make(map[string]bool)
+
+	for _, line := range lines {
+		line = strings.TrimSpace(line)
+		if line == "" {
+			continue
+		}
+
+		matches := lineRe.FindStringSubmatch(line)
+		if matches == nil {
+			continue
+		}
+
+		file := matches[1]
+		lineNum, _ := strconv.Atoi(matches[2])
+		context := strings.TrimSpace(matches[3])
+
+		// Skip definition
+		if strings.Contains(file, "obiseq/attributes.go") && lineNum == 132 {
+			continue
+		}
+
+		// Extract key
+		var key string
+		if keyMatches := keyRe.FindStringSubmatch(context); keyMatches != nil {
+			key = keyMatches[1]
+		} else if tmplMatches := templateKeyRe.FindStringSubmatch(context); tmplMatches != nil {
+			key = tmplMatches[1]
+		} else {
+			continue
+		}
+
+		// Get function name using treesitter
+		funcName := getFunctionNameTreesitter(file, lineNum)
+
+		uniqueKey := fmt.Sprintf("%s:%d", file, lineNum)
+		if seen[uniqueKey] {
+			continue
+		}
+		seen[uniqueKey] = true
+
+		refs = append(refs, Reference{
+			File:     filepath.Base(file),
+			Line:     lineNum,
+			Column:   0,
+			Key:      key,
+			Function: funcName,
+			Context:  context,
+		})
+	}
+
+	sort.Slice(refs, func(i, j int) bool {
+		if refs[i].File != refs[j].File {
+			return refs[i].File < refs[j].File
+		}
+		return refs[i].Line < refs[j].Line
+	})
+
+	result := Result{
+		Method:     "SetAttribute",
+		Signature:  "func (s *BioSequence) SetAttribute(key string, value interface{})",
+		Definition: basePath + "/pkg/obiseq/attributes.go:132",
+		References: refs,
+		Total:      len(refs),
+	}
+
+	outputJSON, err := json.MarshalIndent(result, "", "  ")
+	if err != nil {
+		fmt.Fprintf(os.Stderr, "Error marshaling JSON: %v\n", err)
+		os.Exit(1)
+	}
+
+	fmt.Println(string(outputJSON))
+}
+
+// getFunctionNameTreesitter uses the treesitter_cursor_walk tool to get the containing function
+func getFunctionNameTreesitter(file string, targetLine int) string {
+	// Convert to 0-based for treesitter
+	row := targetLine - 1
+
+	// Use treesitter cursor walk to get ancestors
+	cmd := exec.Command("bash", "-c",
+		fmt.Sprintf(`kilo treesitter_cursor_walk --file_path %q --row %d --column 0 --max_depth 10 2>/dev/null`, file, row))
+
+	output, err := cmd.Output()
+	if err != nil {
+		return findContainingFunction(file, targetLine)
+	}
+
+	// Parse the JSON output to find function_declaration or method_declaration
+	var result map[string]interface{}
+	if err := json.Unmarshal(output, &result); err != nil {
+		return findContainingFunction(file, targetLine)
+	}
+
+	// Check ancestors for function declaration
+	if ancestors, ok := result["ancestors"].([]interface{}); ok {
+		for _, a := range ancestors {
+			if anc, ok := a.(map[string]interface{}); ok {
+				nodeType, _ := anc["type"].(string)
+				if nodeType == "function_declaration" || nodeType == "method_declaration" {
+					// Try to get the function name from children
+					if children, ok := anc["children"].([]interface{}); ok {
+						for _, c := range children {
+							if child, ok := c.(map[string]interface{}); ok {
+								childType, _ := child["type"].(string)
+								if childType == "identifier" {
+									if text, ok := child["text"].(string); ok {
+										return text
+									}
+								}
+								if childType == "field_identifier" {
+									if text, ok := child["text"].(string); ok {
+										return text
+									}
+								}
+							}
+						}
+					}
+				}
+				if nodeType == "func_literal" {
+					return "closure"
+				}
+			}
+		}
+	}
+
+	return findContainingFunction(file, targetLine)
+}
+
+func findContainingFunction(file string, targetLine int) string {
+	data, err := os.ReadFile(file)
+	if err != nil {
+		return ""
+	}
+	lines := strings.Split(string(data), "\n")
+
+	for i := targetLine - 1; i >= 0 && i >= targetLine-200; i-- {
+		if i >= len(lines) {
+			continue
+		}
+		line := strings.TrimSpace(lines[i])
+
+		if line == "}" && i > 0 {
+			for j := i - 1; j >= 0 && j >= i-50; j-- {
+				if j >= len(lines) {
+					continue
+				}
+				funcLine := strings.TrimSpace(lines[j])
+				if strings.HasPrefix(funcLine, "func ") {
+					if match := regexp.MustCompile(`func\s+\([^)]+\)\s+([a-zA-Z_][a-zA-Z0-9_]*)\s*\(`).FindStringSubmatch(funcLine); match != nil {
+						return match[1]
+					}
+					if match := regexp.MustCompile(`func\s+([a-zA-Z_][a-zA-Z0-9_]*)\s*\(`).FindStringSubmatch(funcLine); match != nil {
+						return match[1]
+					}
+				}
+			}
+			continue
+		}
+
+		if strings.HasPrefix(line, "func ") {
+			if match := regexp.MustCompile(`func\s+\([^)]+\)\s+([a-zA-Z_][a-zA-Z0-9_]*)\s*\(`).FindStringSubmatch(line); match != nil {
+				return match[1]
+			}
+			if match := regexp.MustCompile(`func\s+([a-zA-Z_][a-zA-Z0-9_]*)\s*\(`).FindStringSubmatch(line); match != nil {
+				return match[1]
+			}
+		}
+	}
+
+	return ""
+}

scripts/find_setattribute.sh

@@ -0,0 +1,36 @@
+#!/bin/bash
+
+basePath="/Users/coissac/Sync/travail/__MOI__/GO/obitools4"
+OUTPUT_FILE="${1:-/dev/stdout}"
+
+# Get all SetAttribute calls
+rg -n '\.SetAttribute\(' "$basePath/pkg" --type go | while read -r line; do
+	file="${line%%:*}"
+	line_num="${line%:*}"
+	line_num="${line_num##*:}"
+	context="${line##*: }"
+
+	# Extract key (only literal strings)
+	key=$(echo "$context" | sed -n 's/.*SetAttribute("\([^"]*\)".*/\1/p')
+	[ -z "$key" ] && continue
+
+	# Get function name using treesitter
+	func=$(kilo treesitter_cursor_walk \
+		--file_path "$file" \
+		--row "$((line_num - 1))" \
+		--column 0 \
+		--max_depth 10 2>/dev/null |
+		jq -r '.ancestors[] | select(.type == "function_declaration" or .type == "method_declaration") | .children[] | select(.type == "identifier" or .type == "field_identifier") | .text' 2>/dev/null)
+
+	# Fallback to func_literal for closures
+	if [ -z "$func" ]; then
+		func=$(kilo treesitter_cursor_walk \
+			--file_path "$file" \
+			--row "$((line_num - 1))" \
+			--column 0 \
+			--max_depth 10 2>/dev/null |
+			jq -r '.ancestors[] | select(.type == "func_literal") | "closure"' 2>/dev/null)
+	fi
+
+	echo "$(basename "$file")|$line_num|$key|${func:-unknown}|$context"
+done | sort -t'|' -k1,1 -k2,2n

setattribute_refs.json

@@ -0,0 +1,308 @@
+{
+  "obiannotate": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ],
+    "(git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiiter.IBioSequence).NumberSequences$1": [
+      "seq_number"
+    ]
+  },
+  "obiclean": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obicleandb": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetCount": [
+      "count"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obicomplement": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obiconsensus": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetCount": [
+      "count"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ],
+    "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obiconsensus.BuildConsensus": [
+      "obiconsensus_kmer_max_occur",
+      "obiconsensus_filtered_graph_size",
+      "obiconsensus_full_graph_size",
+      "obiconsensus_consensus",
+      "obiconsensus_weight",
+      "obiconsensus_seq_length",
+      "obiconsensus_kmer_size"
+    ],
+    "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obiconsensus.MinionClusterDenoise": [
+      "obiconsensus_consensus"
+    ],
+    "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obiconsensus.MinionDenoise$1": [
+      "obiconsensus_consensus",
+      "obiconsensus_weight"
+    ]
+  },
+  "obiconvert": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obicount": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obicsv": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obidemerge": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obidistribute": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obigrep": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obijoin": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obikmermatch": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obikmersimcount": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obilandmark": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetCoordinate": [
+      "landmark_coord"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetOBITagGeomRefIndex": [
+      "obitag_geomref_index"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ],
+    "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obilandmark.CLISelectLandmarkSequences": [
+      "landmark_id"
+    ]
+  },
+  "obimatrix": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obimicrosat": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obimultiplex": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obipairing": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence)._revcmpMutation": [
+      "pairing_mismatches"
+    ],
+    "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obialign.BuildQualityConsensus": [
+      "pairing_mismatches"
+    ]
+  },
+  "obipcr": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obireffamidx": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetOBITagRefIndex": [
+      "obitag_ref_index"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ],
+    "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitools/obirefidx.IndexFamilyDB": [
+      "reffamidx_id"
+    ]
+  },
+  "obirefidx": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetOBITagRefIndex": [
+      "obitag_ref_index"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obiscript": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obisplit": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obisummary": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obitag": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetPath": [
+      "taxonomic_path"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  },
+  "obitagpcr": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obingslibrary.NGSLibrary).ExtractMultiBarcode": [
+      "obimultiplex_error",
+      "obimultiplex_amplicon_rank"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence)._revcmpMutation": [
+      "pairing_mismatches"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence)._subseqMutation": [
+      "pairing_mismatches"
+    ],
+    "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obialign.BuildQualityConsensus": [
+      "pairing_mismatches"
+    ]
+  },
+  "obitaxonomy": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetPath": [
+      "taxonomic_path"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ],
+    "(git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiiter.IBioSequence).NumberSequences$1": [
+      "seq_number"
+    ]
+  },
+  "obiuniq": {
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetCount": [
+      "count"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetDefinition": [
+      "definition"
+    ],
+    "(*git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq.BioSequence).SetTaxid": [
+      "taxid"
+    ]
+  }
+}

version.txt

@@ -1 +1 @@
-4.4.24
+4.5.0

x

@@ -0,0 +1,19 @@
+```markdown
+# DNA Scoring and Matching Utilities in `obialign`
+
+This module provides low-level utilities for computing nucleotide alignment scores using probabilistic and bit-encoded representations.
+
+- **Bit Encoding**: Nucleotides are encoded in 4-bit groups (e.g., `A=0b0001`, `C=0b0010`, etc.), enabling efficient bitwise comparison.
+- **`_MatchRatio(a, b)`**: Computes a normalized match ratio between two encoded bytes based on shared bits:
+  `ratio = common_bits / (bits_in_a × bits_in_b)`.
+- **`_FourBitsCount`**: Precomputed lookup table for Hamming weight (popcount) of 4-bit values.
+- **Log-space Arithmetic**: Helper functions (`_Logaddexp`, `_Logdiffexp`, `_Log1mexp`) ensure numerical stability in probabilistic computations.
+- **Phred-scaled Quality Integration**:
+  `_MatchScoreRatio(QF, QR)` derives log-odds match/mismatch scores from Phred quality values (`QF`, `QR`), modeling sequencing error probabilities.
+- **Precomputed Matrices**:
+  - `_NucPartMatch[i][j]`: Match ratios for all nucleotide pairs (from 4-bit codes).
+  - `_NucScorePartMatchMatch/Mismatch[i][j]`: Integer-scaled match/mismatch scores (×10) for quality pairs `(i, j)` in `[0..99]`.
+- **Thread-Safe Initialization**: `_InitDNAScoreMatrix()` ensures one-time, synchronized initialization of all scoring tables via a mutex.
+
+Designed for high-performance alignment kernels where speed and numerical robustness are critical.
+```