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.

.gitignore

@@ -21,6 +21,8 @@ xx
 
 .rhistory
 /.vscode
+
+/benchmarck
 /build
 /bugs
 autodoc

.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 {

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,16 +496,26 @@ func (pattern ApatPattern) AllMatches(sequence ApatSequence, begin, length int)
 				(*cpattern)[0:int(pattern.pointer.pointer.patlen)],
 				frg)
 
-			// olderr := m[2]
+			if pb < 0 {
-			m[2] = score
+				// obialign.LocatePattern could not reconstruct a
-			m[0] = start + pb
+				// reliable alignment (e.g. the match sits too close
-			m[1] = start + pe
+				// 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
+				m[1] = start + pe
 
-			// obilog.Warnf("seq[%d@%d:%d] %d: %s %d - %s:%s:%s", i, m[0], m[1], olderr, sequence.pointer.reference.Id(), score,
+				// 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]])
+				// 	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,8 +199,111 @@ 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()
+	var err error
+
+	skey := obiutils.UnsafeString(key)
+
+	switch {
+	case skey == "id":
+		sequence.SetId(string(value))
+	case skey == "definition":
+		sequence.SetDefinition(string(value))
+
+	case skey == "count":
+		if dataType != jsonparser.Number {
+			log.Fatalf("%s: Count attribut must be numeric: %s", sequence.Id(), string(value))
+		}
+		count, err := jsonparser.ParseInt(value)
+		if err != nil {
+			log.Fatalf("%s: Cannot parse count %s", sequence.Id(), string(value))
+		}
+		sequence.SetCount(int(count))
+
+	case skey == "obiclean_weight":
+		weight, err := _parse_json_map_int(value)
+		if err != nil {
+			log.Fatalf("%s: Cannot parse obiclean weight %s", sequence.Id(), string(value))
+		}
+		annotations[skey] = weight
+
+	case skey == "obiclean_status":
+		status, err := _parse_json_map_string(value)
+		if err != nil {
+			log.Fatalf("%s: Cannot parse obiclean status %s", sequence.Id(), string(value))
+		}
+		annotations[skey] = status
+
+	case strings.HasPrefix(skey, "merged_"):
+		if dataType == jsonparser.Object {
+			data, err := _parse_json_map_int(value)
+			if err != nil {
+				log.Fatalf("%s: Cannot parse merged slot %s: %v", sequence.Id(), skey, err)
+			} else {
+				annotations[skey] = obiseq.MapAsStatsOnValues(data)
+			}
+		} else {
+			log.Fatalf("%s: Cannot parse merged slot %s", sequence.Id(), skey)
+		}
+
+	case skey == "taxid":
+		if dataType == jsonparser.Number || dataType == jsonparser.String {
+			taxid := string(value)
+			sequence.SetTaxid(taxid)
+		} else {
+			log.Fatalf("%s: Cannot parse taxid %s", sequence.Id(), string(value))
+		}
+
+	case strings.HasSuffix(skey, "_taxid"):
+		if dataType == jsonparser.Number || dataType == jsonparser.String {
+			rank := skey[:len(skey)-len("_taxid")]
+
+			taxid := string(value)
+			sequence.SetTaxid(taxid, rank)
+		} else {
+			log.Fatalf("%s: Cannot parse taxid %s", sequence.Id(), string(value))
+		}
+
+	default:
+		skey = strings.Clone(skey)
+		switch dataType {
+		case jsonparser.String:
+			annotations[skey] = string(value)
+		case jsonparser.Number:
+			// Try to parse the number as an int at first then as float if that fails.
+			annotations[skey], err = jsonparser.ParseInt(value)
+			if err != nil {
+				annotations[skey], err = strconv.ParseFloat(obiutils.UnsafeString(value), 64)
+			}
+		case jsonparser.Array:
+			annotations[skey], err = _parse_json_array_interface(value)
+		case jsonparser.Object:
+			annotations[skey], err = _parse_json_map_interface(value)
+		case jsonparser.Boolean:
+			annotations[skey], err = jsonparser.ParseBoolean(value)
+		case jsonparser.Null:
+			annotations[skey] = nil
+		default:
+			log.Fatalf("Unknown data type %v", dataType)
+		}
+	}
+
+	if err != nil {
+		annotations[skey] = "NaN"
+		log.Fatalf("%s: Cannot parse value %s assicated to key %s into a %s value",
+			sequence.Id(), string(value), skey, dataType.String())
+	}
+
+	return err
+}
+
+func _parse_json_header_(header string, sequence *obiseq.BioSequence) string {
 	start := -1
 	stop := -1
 	level := 0
@@ -240,101 +343,7 @@ func _parse_json_header_(header string, sequence *obiseq.BioSequence) string {
 
 	jsonparser.ObjectEach(obiutils.UnsafeBytes(header[start:stop]),
 		func(key []byte, value []byte, dataType jsonparser.ValueType, offset int) error {
-			var err error
+			return _parse_json_annotation_field(key, value, dataType, sequence)
-
-			skey := obiutils.UnsafeString(key)
-
-			switch {
-			case skey == "id":
-				sequence.SetId(string(value))
-			case skey == "definition":
-				sequence.SetDefinition(string(value))
-
-			case skey == "count":
-				if dataType != jsonparser.Number {
-					log.Fatalf("%s: Count attribut must be numeric: %s", sequence.Id(), string(value))
-				}
-				count, err := jsonparser.ParseInt(value)
-				if err != nil {
-					log.Fatalf("%s: Cannot parse count %s", sequence.Id(), string(value))
-				}
-				sequence.SetCount(int(count))
-
-			case skey == "obiclean_weight":
-				weight, err := _parse_json_map_int(value)
-				if err != nil {
-					log.Fatalf("%s: Cannot parse obiclean weight %s", sequence.Id(), string(value))
-				}
-				annotations[skey] = weight
-
-			case skey == "obiclean_status":
-				status, err := _parse_json_map_string(value)
-				if err != nil {
-					log.Fatalf("%s: Cannot parse obiclean status %s", sequence.Id(), string(value))
-				}
-				annotations[skey] = status
-
-			case strings.HasPrefix(skey, "merged_"):
-				if dataType == jsonparser.Object {
-					data, err := _parse_json_map_int(value)
-					if err != nil {
-						log.Fatalf("%s: Cannot parse merged slot %s: %v", sequence.Id(), skey, err)
-					} else {
-						annotations[skey] = obiseq.MapAsStatsOnValues(data)
-					}
-				} else {
-					log.Fatalf("%s: Cannot parse merged slot %s", sequence.Id(), skey)
-				}
-
-			case skey == "taxid":
-				if dataType == jsonparser.Number || dataType == jsonparser.String {
-					taxid := string(value)
-					sequence.SetTaxid(taxid)
-				} else {
-					log.Fatalf("%s: Cannot parse taxid %s", sequence.Id(), string(value))
-				}
-
-			case strings.HasSuffix(skey, "_taxid"):
-				if dataType == jsonparser.Number || dataType == jsonparser.String {
-					rank := skey[:len(skey)-len("_taxid")]
-
-					taxid := string(value)
-					sequence.SetTaxid(taxid, rank)
-				} else {
-					log.Fatalf("%s: Cannot parse taxid %s", sequence.Id(), string(value))
-				}
-
-			default:
-				skey = strings.Clone(skey)
-				switch dataType {
-				case jsonparser.String:
-					annotations[skey] = string(value)
-				case jsonparser.Number:
-					// Try to parse the number as an int at first then as float if that fails.
-					annotations[skey], err = jsonparser.ParseInt(value)
-					if err != nil {
-						annotations[skey], err = strconv.ParseFloat(obiutils.UnsafeString(value), 64)
-					}
-				case jsonparser.Array:
-					annotations[skey], err = _parse_json_array_interface(value)
-				case jsonparser.Object:
-					annotations[skey], err = _parse_json_map_interface(value)
-				case jsonparser.Boolean:
-					annotations[skey], err = jsonparser.ParseBoolean(value)
-				case jsonparser.Null:
-					annotations[skey] = nil
-				default:
-					log.Fatalf("Unknown data type %v", dataType)
-				}
-			}
-
-			if err != nil {
-				annotations[skey] = "NaN"
-				log.Fatalf("%s: Cannot parse value %s assicated to key %s into a %s value",
-					sequence.Id(), string(value), skey, dataType.String())
-			}
-
-			return err
 		},
 	)
 

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,28 +301,44 @@ func ParseOBIFeatures(text string, annotations obiseq.Annotation) string {
 					stop = m[1] + 1
 				} else {
 
-					// Generic value
+					// array value
-
+					m = __match__array__(part)
-					// m = __obi_header_value_general_pattern__.FindIndex(part)
-					m = __match__general__(part)
 					if len(m) > 0 {
 						bvalue = bytes.TrimSpace(part[m[0]:(m[1] - 1)])
-
+						j := bytes.ReplaceAll(bvalue, []byte("'"), []byte(`"`))
-						if __is_false__(bvalue) {
+						j = __obi_header_map_int_key__.ReplaceAll(j, []byte(`$1"$2":`))
-							value = false
+						arr, err := _parse_json_array_interface(j)
+						if err != nil {
+							value = string(bvalue)
 						} else {
-							if __is_true__(bvalue) {
+							value = arr
-								value = true
-							} else {
-								value = string(bvalue)
-							}
 						}
-
 						stop = m[1] + 1
 					} else {
-						// no value
+
-						break
+						// Generic value
-					} // End of No value
+
+						// m = __obi_header_value_general_pattern__.FindIndex(part)
+						m = __match__general__(part)
+						if len(m) > 0 {
+							bvalue = bytes.TrimSpace(part[m[0]:(m[1] - 1)])
+
+							if __is_false__(bvalue) {
+								value = false
+							} else {
+								if __is_true__(bvalue) {
+									value = true
+								} else {
+									value = string(bvalue)
+								}
+							}
+
+							stop = m[1] + 1
+						} else {
+							// 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,19 +404,19 @@ func WriteFastSeqOBIHeade(buffer *bytes.Buffer, sequence *obiseq.BioSequence) {
 					buffer.WriteString(fmt.Sprintf("%s=", key))
 					buffer.Write(tv)
 					buffer.WriteString("; ")
-				case map[string]int,
-					map[string]string,
-					map[string]interface{}:
-					tv, err := obiutils.JsonMarshal(t)
-					if err != nil {
-						log.Fatalf("Cannot convert %v value", value)
-					}
-					tv = bytes.ReplaceAll(tv, []byte(`"`), []byte("'"))
-					buffer.WriteString(fmt.Sprintf("%s=", key))
-					buffer.Write(tv)
-					buffer.WriteString("; ")
 				default:
-					buffer.WriteString(fmt.Sprintf("%s=%v; ", key, value))
+					if obiutils.IsAMap(value) || obiutils.IsASlice(value) || obiutils.IsAnArray(value) {
+						tv, err := obiutils.JsonMarshal(t)
+						if err != nil {
+							log.Fatalf("Cannot convert %v value", value)
+						}
+						tv = bytes.ReplaceAll(tv, []byte(`"`), []byte("'"))
+						buffer.WriteString(fmt.Sprintf("%s=", key))
+						buffer.Write(tv)
+						buffer.WriteString("; ")
+					} else {
+						buffer.WriteString(fmt.Sprintf("%s=%v; ", key, value))
+					}
 				}
 			}
 		}

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/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
 		}
@@ -129,24 +147,20 @@ func KmerEntropy(kmer uint64, k int, levelMax int) float64 {
 // IMPORTANT: a KmerEntropyFilter is NOT safe for concurrent use.
 // Each goroutine must create its own instance via NewKmerEntropyFilter.
 type KmerEntropyFilter struct {
-	k          int
+	k                   int
-	levelMax   int
+	levelMax            int
-	threshold  float64
+	threshold           float64
-	nLogN      []float64
+	nLogN               []float64
-	normTables [][]int
+	normTables          [][]int
-	emaxValues []float64
+	classLogSizeTables  [][]float64
-	logNwords  []float64
+	emaxValues          []float64
+	logNwords           []float64
 	// Pre-allocated frequency tables reused across Entropy() calls.
 	// One per word size (index 0 unused). Reset to zero before each use.
 	freqTables [][]int
 }
 
 // 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))
 		}
@@ -195,14 +227,15 @@ func NewKmerEntropyFilter(k, levelMax int, threshold float64) *KmerEntropyFilter
 	}
 
 	return &KmerEntropyFilter{
-		k:          k,
+		k:                  k,
-		levelMax:   levelMax,
+		levelMax:           levelMax,
-		threshold:  threshold,
+		threshold:          threshold,
-		nLogN:      nLogN,
+		nLogN:              nLogN,
-		normTables: normTables,
+		normTables:         normTables,
-		emaxValues: emaxValues,
+		classLogSizeTables: classLogSizeTables,
-		logNwords:  logNwords,
+		emaxValues:         emaxValues,
-		freqTables: freqTables,
+		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/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

@@ -6,13 +6,14 @@ import (
 )
 
 type _Options struct {
-	discardErrors   bool
+	discardErrors      bool
-	unidentified    string
+	unidentified       string
-	allowedMismatch int
+	allowedMismatch    int
-	allowsIndel     bool
+	allowedMismatchSet bool
-	withProgressBar bool
+	allowsIndel        bool
-	parallelWorkers int
+	withProgressBar    bool
-	batchSize       int
+	parallelWorkers    int
+	batchSize          int
 }
 
 // Options stores a set of option usable by the
@@ -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.40"
+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/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/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 {

version.txt

@@ -1 +1 @@
-4.4.40
+4.5.0