Add benchmark pipeline for dense and sparse query testing
Introduces a complete query benchmark track to evaluate performance and verify consistency between dense and sparse index formats. Adds scripts to simulate fixed-size paired-end reads, pack a sparse presence index, execute queries in both modes, and capture wall time and RSS metrics. Includes a verification step that compares outputs by read ID to ensure content identity across parallel processing. Updates build configuration, documentation, and ignore patterns to support the new pipeline for two microbial specimens.
This commit is contained in:
Eric Coissac
+12
-1
@@ -10,13 +10,23 @@ src/target
|
||||
data-stress
|
||||
*.fasta
|
||||
*.fasta.gz
|
||||
*.fastq
|
||||
*.fastq.gz
|
||||
*.vcf
|
||||
*.zst
|
||||
*.zst.meta
|
||||
*.pb
|
||||
./**/*.json
|
||||
**/*.json
|
||||
*.bin
|
||||
*.log
|
||||
*.csv
|
||||
*.meta
|
||||
*.pfiv
|
||||
*.done
|
||||
*.efh
|
||||
*.efl
|
||||
*.bin.idx
|
||||
*.prsb
|
||||
Betula_exilis--IGA-24-33
|
||||
benchmark/genomes
|
||||
benchmark/genomes_orig
|
||||
@@ -34,6 +44,7 @@ benchmark/reference_dist
|
||||
benchmark/obikmer_dist
|
||||
benchmark/specific_index_count
|
||||
benchmark/specific_index_presence
|
||||
benchmark/query_data
|
||||
TNT
|
||||
phyg
|
||||
biblio
|
||||
|
||||
@@ -1,84 +1,93 @@
|
||||
# Benchmark: query-path testing (discussion)
|
||||
# Benchmark: query-path testing
|
||||
|
||||
`benchmark/Makefile` exercises indexing, merge, and phylo distance
|
||||
reconstruction against simulated bacterial genomes, but has no coverage of
|
||||
`obikmer query` — the read-matching path — nor of the sparse packed
|
||||
presence-matrix format (`obikmer pack --sparse`). This note captures the
|
||||
planned extension.
|
||||
reconstruction against simulated bacterial genomes. It now also covers
|
||||
`obikmer query` — the read-matching path — and the sparse packed
|
||||
presence-matrix format (`obikmer pack --sparse`), previously untested by
|
||||
this pipeline.
|
||||
|
||||
## Motivation
|
||||
|
||||
- `query` is untested end-to-end. A regression there would not be caught by
|
||||
the existing `verify_presence`/`verify_merge_presence` branches, which only
|
||||
check index *content* against the `.npz` truth, never the query API.
|
||||
- `query` had no end-to-end coverage. A regression there would not be caught
|
||||
by `verify_presence`/`verify_merge_presence`, which only check index
|
||||
*content* against the `.npz` truth, never the query API.
|
||||
- `pack --sparse` produces a presence-matrix format documented (see
|
||||
[siblings.md](../architecture/siblings.md)) as faster for single-row
|
||||
access (query) and slower for column-oriented access (phylo `--metric`).
|
||||
`global_index_presence/` built by `merge_presence.sh` is currently always
|
||||
packed dense (packing is a stage inside `merge`, not a separate `pack`
|
||||
invocation). There is no dense/sparse regression check.
|
||||
`global_index_presence/` built by `merge_presence.sh` is always packed
|
||||
dense (packing is a stage inside `merge`, not a separate `pack`
|
||||
invocation) — there was no dense/sparse regression check.
|
||||
|
||||
## Plan
|
||||
## Query read source
|
||||
|
||||
**New read source, independent of `simulated_data/`.** Reusing
|
||||
`simulated_data/<species>/<strain>/reads_R1.fastq.gz` for queries would bias
|
||||
the test: those reads were already folded into the index being queried, with
|
||||
the same sequencing-error draw. Query reads must come from a *second*,
|
||||
independent `iss generate` run against the same reference genome(s) — new
|
||||
random error draw, same underlying sequence — landing in a separate tree:
|
||||
`query_data/<species>/<strain>/reads_R1.fastq.gz`, built by the existing
|
||||
`simulate_one.sh` (unseeded, so a second invocation naturally draws different
|
||||
reads).
|
||||
Query reads are independent of `simulated_data/` (which is folded into the
|
||||
index being queried): reusing those reads would test against the exact
|
||||
error draw the index was built from. `query_data/<species>/<strain>/` holds
|
||||
a *second*, independent `iss generate` run against the same reference
|
||||
genome, via `simulate_query_one.sh` — unseeded, so a second draw picks up
|
||||
different sequencing errors than `simulate_one.sh`'s draw for the same
|
||||
genome. Fixed at 100,000 read pairs per genome (not coverage-proportional
|
||||
like the 15x used for `simulated_data/`), so wall/RSS numbers stay
|
||||
comparable across genomes of very different sizes.
|
||||
|
||||
Two specimens chosen as query sources (enough to catch a dense/sparse
|
||||
regression without duplicating the exhaustive per-specimen coverage
|
||||
`verify_merge_presence` already provides across all `SPECIMENS`):
|
||||
`Escherichia_coli--K-12_MG1655` (common, well-represented bacterium) and
|
||||
`Saccharolobus_islandicus--M.16.4` (the only archaeon in `SPECIES` — distant
|
||||
lineage, different GC content, stresses the query path differently from a
|
||||
close-relative match).
|
||||
Two query-source specimens, hardcoded as `QUERY_SPECIMENS` in
|
||||
`make_deps.py`: `Escherichia_coli--K-12_MG1655` (common, well-represented
|
||||
bacterium) and `Saccharolobus_islandicus--M.16.4` (the only archaeon in
|
||||
`SPECIES` — distant lineage, stresses the query path differently from a
|
||||
close-relative match). Two is enough to catch a dense/sparse regression
|
||||
without duplicating the exhaustive per-specimen coverage
|
||||
`verify_merge_presence` already provides across all `SPECIMENS`.
|
||||
|
||||
`make_deps.py` needs a `QUERY_SPECIMENS` list (explicit, short) and, for each,
|
||||
an extra dependency line:
|
||||
```
|
||||
query_data/<species>/<strain>/reads_R1.fastq.gz: genomes/<genome>.fna.gz
|
||||
```
|
||||
distinct from the `simulated_data/...` rule for the same specimen.
|
||||
## Sparse global index
|
||||
|
||||
Read count fixed at 100,000 read pairs per genome, independent of genome
|
||||
size — unlike `simulate_one.sh`'s `simulated_data/` runs, which derive
|
||||
`n_reads` from a fixed 15x coverage target. A query benchmark does not need
|
||||
coverage-proportional depth; a fixed pair count keeps the two query runs
|
||||
comparable to each other and keeps wall/RSS numbers meaningful across
|
||||
genomes of very different sizes (bacterium vs archaeon). This likely needs a
|
||||
dedicated `simulate_query_one.sh` (or a parameter to `simulate_one.sh`)
|
||||
rather than reusing it unchanged, since `n_reads` is currently computed
|
||||
in-script from genome size.
|
||||
`global_index_presence_sparse/` is built by `pack_sparse.sh`: copy
|
||||
`global_index_presence/` wholesale, then `obikmer pack --sparse` in place.
|
||||
This works directly because `merge`'s pack stage (`merge.rs:252`,
|
||||
`pack_matrices(false)`) keeps the per-genome column files on disk after
|
||||
dense-packing — `pack_sparse_bit_matrix` (`obicompactvec/src/bitmatrix/sparse.rs:447`)
|
||||
reads those, is idempotent, and removes `matrix.pbmx` once the sparse form
|
||||
is written, so `Persistent::open` falls through to the sparse format
|
||||
afterward. No separate merge run needed.
|
||||
|
||||
**Phase 1 — sparse global index.** New target
|
||||
`global_index_presence_sparse/index.done`, built from `global_index_presence/`
|
||||
via `obikmer pack --sparse`. Open question, to verify against the `pack`
|
||||
implementation before writing the rule: does `pack --sparse` accept an
|
||||
already dense-packed index in place (`cp -r` + repack), or does it require
|
||||
the pre-pack column layout, forcing a dedicated merge run instead of reusing
|
||||
`global_index_presence/`?
|
||||
## Query runs
|
||||
|
||||
**Phase 2 — query runs.** For each of the two `QUERY_SPECIMENS`, run
|
||||
`obikmer query` against both `global_index_presence` and
|
||||
`global_index_presence_sparse`, capturing Reporter wall/RSS stats the same
|
||||
way `merge_presence.sh` does (stderr capture + `parse_reporter`).
|
||||
`query_one.sh dense|sparse SPECIMEN` runs `obikmer query --count-missing`
|
||||
against `global_index_presence` or `global_index_presence_sparse`, output
|
||||
gzipped to `query_{dense,sparse}/SPECIMEN.fasta.gz`, Reporter wall/RSS
|
||||
captured to `stats/query_{dense,sparse}/SPECIMEN.stats` (same
|
||||
stderr-parsing convention as `merge_presence.sh`).
|
||||
|
||||
**Phase 3 — dense/sparse regression.** `verify_query.py` diffs the two query
|
||||
JSON outputs per specimen (same matches, same per-genome presence
|
||||
annotations) → `.stats` CSV (`run,specimen,mismatches,pct`), aggregated by
|
||||
`aggregate_stats.sh` under a new `query` case. Any mismatch is a real
|
||||
regression — dense and sparse must be content-identical, only I/O access
|
||||
pattern differs.
|
||||
Flags: `--count-missing` only. `--mismatch` is a no-op today
|
||||
(`query/mod.rs:212-213`, prints "not yet implemented, ignored") — left off
|
||||
rather than tested for a feature that doesn't exist yet.
|
||||
|
||||
**Phase 4 — performance comparison.** No dedicated script: the wall/RSS
|
||||
columns from Phase 2's `.stats` files, aggregated, are the dense-vs-sparse
|
||||
performance comparison (the expected win for query on sparse, per the `pack
|
||||
--sparse` help text).
|
||||
## Dense/sparse regression
|
||||
|
||||
`count` track is out of scope for the sparse branch: `pack --sparse` targets
|
||||
presence matrices only (per CLI help), no count equivalent confirmed.
|
||||
`verify_query.py` compares the two query outputs per specimen, matched by
|
||||
read id (not stream position — the query pipeline chunks input across
|
||||
worker threads and doesn't guarantee output order). Compares `kmer_count`,
|
||||
`kmer_missing`, and the full `kmer_strict_matches` map per read. Any
|
||||
mismatch is a real regression: dense and sparse must be content-identical,
|
||||
only I/O access pattern differs. `.stats` → `stats/verify_query/`,
|
||||
aggregated by `aggregate_stats.sh query|verify_query`-style cases
|
||||
(`query_dense`, `query_sparse`, `verify_query`).
|
||||
|
||||
## Performance comparison
|
||||
|
||||
No dedicated script: the wall/RSS columns from the `query_dense` and
|
||||
`query_sparse` aggregated `.stats` CSVs are the dense-vs-sparse performance
|
||||
comparison — the expected win for query on sparse, per the `pack --sparse`
|
||||
help text.
|
||||
|
||||
## Scope
|
||||
|
||||
`count` track excluded from the sparse branch: `pack --sparse` targets
|
||||
presence matrices only (per CLI help); `pack_matrices` leaves count
|
||||
matrices untouched regardless of the `sparse` flag
|
||||
(`obikindex/src/index.rs:308`).
|
||||
|
||||
## New Makefile targets
|
||||
|
||||
`simulate_query`, `pack_sparse`, `query_dense`, `query_sparse`,
|
||||
`aggregate_query_dense`, `aggregate_query_sparse`, `verify_query`,
|
||||
`aggregate_verify_query` — the last three folded into `all`.
|
||||
|
||||
+60
-4
@@ -38,6 +38,12 @@ SPECIFIC_PRESENCE_STATS := $(SPECIES:%=stats/specific_kmer_presence/%.stats)
|
||||
SPECIFIC_COUNT_DONE := $(SPECIES:%=specific_index_count/%/index.done)
|
||||
SPECIFIC_COUNT_STATS := $(SPECIES:%=stats/specific_kmer_count/%.stats)
|
||||
SIMULATED_READS := $(foreach s,$(SPECIMENS),simulated_data/$(subst --,/,$s)/reads_R1.fastq.gz)
|
||||
QUERY_READS := $(foreach s,$(QUERY_SPECIMENS),query_data/$(subst --,/,$s)/reads_R1.fastq.gz)
|
||||
QUERY_DENSE_DONE := $(QUERY_SPECIMENS:%=query_dense/%.fasta.gz)
|
||||
QUERY_DENSE_STATS := $(QUERY_SPECIMENS:%=stats/query_dense/%.stats)
|
||||
QUERY_SPARSE_DONE := $(QUERY_SPECIMENS:%=query_sparse/%.fasta.gz)
|
||||
QUERY_SPARSE_STATS := $(QUERY_SPECIMENS:%=stats/query_sparse/%.stats)
|
||||
VERIFY_QUERY_STATS := $(QUERY_SPECIMENS:%=stats/verify_query/%.stats)
|
||||
|
||||
.NOTPARALLEL:
|
||||
|
||||
@@ -51,7 +57,11 @@ SIMULATED_READS := $(foreach s,$(SPECIMENS),simulated_data/$(subst --,/,$s)/read
|
||||
aggregate_verify_presence aggregate_verify_count \
|
||||
verify_merge_presence verify_merge_count \
|
||||
filter_presence filter_count \
|
||||
aggregate_filter_presence aggregate_filter_count
|
||||
aggregate_filter_presence aggregate_filter_count \
|
||||
pack_sparse simulate_query \
|
||||
query_dense query_sparse \
|
||||
aggregate_query_dense aggregate_query_sparse \
|
||||
verify_query aggregate_verify_query
|
||||
|
||||
verify_merge_presence: stats/verify_merge_presence/current.csv
|
||||
verify_merge_count: stats/verify_merge_count/current.csv
|
||||
@@ -59,12 +69,13 @@ verify_merge_count: stats/verify_merge_count/current.csv
|
||||
all: aggregate_verify_presence aggregate_verify_count \
|
||||
verify_merge_presence verify_merge_count \
|
||||
aggregate_filter_presence aggregate_filter_count \
|
||||
dist_comparison
|
||||
dist_comparison \
|
||||
aggregate_query_dense aggregate_query_sparse aggregate_verify_query
|
||||
|
||||
# ── dependency file ───────────────────────────────────────────────────────────
|
||||
|
||||
deps.mk: $(GENOMES)
|
||||
$(VENV_PY) make_deps.py $^ > $@
|
||||
deps.mk: $(GENOMES) make_deps.py
|
||||
$(VENV_PY) make_deps.py $(GENOMES) > $@
|
||||
|
||||
# ── simulation ────────────────────────────────────────────────────────────────
|
||||
# Prerequisites (genome → reads) are in deps.mk; $< is the genome file.
|
||||
@@ -74,6 +85,14 @@ $(SIMULATED_READS):
|
||||
|
||||
simulate: $(SIMULATED_READS)
|
||||
|
||||
# ── query read simulation (fixed size, independent draw) ───────────────────────
|
||||
# Prerequisites (genome → reads) are in deps.mk; $< is the genome file.
|
||||
|
||||
$(QUERY_READS):
|
||||
bash simulate_query_one.sh $< $(dir $@)
|
||||
|
||||
simulate_query: $(QUERY_READS)
|
||||
|
||||
# ── reference kmer sets ───────────────────────────────────────────────────────
|
||||
# Prerequisites (reads → npz) are in deps.mk.
|
||||
|
||||
@@ -228,3 +247,40 @@ stats/verify_merge_presence/current.csv: $(REF_NPZS) global_index_presence/index
|
||||
|
||||
stats/verify_merge_count/current.csv: $(REF_NPZS) global_index_count/index.done
|
||||
bash verify_merge_count.sh
|
||||
|
||||
# ── sparse presence index (query benchmark) ─────────────────────────────────────
|
||||
|
||||
global_index_presence_sparse/index.done: global_index_presence/index.done $(BINARY)
|
||||
bash pack_sparse.sh
|
||||
|
||||
pack_sparse: global_index_presence_sparse/index.done
|
||||
|
||||
# ── query: dense vs sparse ───────────────────────────────────────────────────────
|
||||
# Prerequisites (reads + index → output + .stats) are in deps.mk.
|
||||
|
||||
query_dense/%.fasta.gz \
|
||||
stats/query_dense/%.stats &: $(BINARY)
|
||||
bash query_one.sh dense $*
|
||||
|
||||
query_sparse/%.fasta.gz \
|
||||
stats/query_sparse/%.stats &: $(BINARY)
|
||||
bash query_one.sh sparse $*
|
||||
|
||||
query_dense: $(QUERY_DENSE_DONE)
|
||||
query_sparse: $(QUERY_SPARSE_DONE)
|
||||
|
||||
aggregate_query_dense: $(QUERY_DENSE_STATS)
|
||||
bash aggregate_stats.sh query_dense
|
||||
|
||||
aggregate_query_sparse: $(QUERY_SPARSE_STATS)
|
||||
bash aggregate_stats.sh query_sparse
|
||||
|
||||
# ── query: dense/sparse regression ──────────────────────────────────────────────
|
||||
|
||||
stats/verify_query/%.stats:
|
||||
bash verify_query_one.sh $*
|
||||
|
||||
verify_query: $(VERIFY_QUERY_STATS)
|
||||
|
||||
aggregate_verify_query: $(VERIFY_QUERY_STATS)
|
||||
bash aggregate_stats.sh verify_query
|
||||
|
||||
+50
-2
@@ -75,10 +75,41 @@ flowchart TD
|
||||
verify_merge_count --> vmc[("stats/verify_merge_count/")]
|
||||
end
|
||||
|
||||
subgraph query ["Query track (2 specimens: E. coli + archaeon)"]
|
||||
GENOMES --> simulate_query
|
||||
simulate_query --> qdata[("query_data/")]
|
||||
|
||||
gpres --> pack_sparse
|
||||
BIN --> pack_sparse
|
||||
pack_sparse --> gsparse[("global_index_presence_sparse/")]
|
||||
|
||||
qdata --> query_dense
|
||||
gpres --> query_dense
|
||||
BIN --> query_dense
|
||||
query_dense --> qd[("query_dense/")]
|
||||
query_dense --> qd_stats[("stats/query_dense/")]
|
||||
qd_stats --> aggregate_query_dense
|
||||
|
||||
qdata --> query_sparse
|
||||
gsparse --> query_sparse
|
||||
BIN --> query_sparse
|
||||
query_sparse --> qs[("query_sparse/")]
|
||||
query_sparse --> qs_stats[("stats/query_sparse/")]
|
||||
qs_stats --> aggregate_query_sparse
|
||||
|
||||
qd --> verify_query
|
||||
qs --> verify_query
|
||||
verify_query --> vq_stats[("stats/verify_query/")]
|
||||
vq_stats --> aggregate_verify_query
|
||||
end
|
||||
|
||||
aggregate_verify_presence --> all
|
||||
aggregate_verify_count --> all
|
||||
vmp --> all
|
||||
vmc --> all
|
||||
aggregate_query_dense --> all
|
||||
aggregate_query_sparse --> all
|
||||
aggregate_verify_query --> all
|
||||
all -. "$(MAKE) re-eval" .-> aggregate_filter_presence
|
||||
all -. "$(MAKE) re-eval" .-> aggregate_filter_count
|
||||
```
|
||||
@@ -105,6 +136,14 @@ flowchart TD
|
||||
| `aggregate_filter_count` | `aggregate_stats.sh` | Aggregate species-specific kmer stats (count) |
|
||||
| `verify_merge_presence` | `verify_merge_presence.sh` | Verify global presence index against all reference sets |
|
||||
| `verify_merge_count` | `verify_merge_count.sh` | Verify global count index against all reference sets |
|
||||
| `simulate_query` | `simulate_query_one.sh` | Simulate a fixed-size (100k pairs) read set per query specimen |
|
||||
| `pack_sparse` | `pack_sparse.sh` | Build `global_index_presence_sparse/` from `global_index_presence/` |
|
||||
| `query_dense` | `query_one.sh dense` | Query each query specimen's reads against the dense global index |
|
||||
| `query_sparse` | `query_one.sh sparse` | Query each query specimen's reads against the sparse global index |
|
||||
| `aggregate_query_dense` | `aggregate_stats.sh` | Aggregate dense query wall/RSS stats |
|
||||
| `aggregate_query_sparse` | `aggregate_stats.sh` | Aggregate sparse query wall/RSS stats |
|
||||
| `verify_query` | `verify_query_one.sh` | Diff dense vs sparse query output per specimen (regression check) |
|
||||
| `aggregate_verify_query` | `aggregate_stats.sh` | Aggregate dense/sparse query regression stats |
|
||||
|
||||
## Directory layout
|
||||
|
||||
@@ -113,13 +152,18 @@ benchmark/
|
||||
├── genomes/ # input reference genomes (.fna.gz)
|
||||
├── simulated_data/ # generated by simulate
|
||||
│ └── <species>/<specimen>/
|
||||
├── query_data/ # generated by simulate_query (2 specimens, fixed 100k pairs)
|
||||
│ └── <species>/<specimen>/
|
||||
├── reference_index/ # reference kmer sets (.npz)
|
||||
├── specimen_index_presence/ # per-specimen presence indexes
|
||||
├── specimen_index_count/ # per-specimen count indexes
|
||||
├── global_index_presence/ # merged global presence index
|
||||
├── global_index_presence/ # merged global presence index (dense-packed)
|
||||
├── global_index_presence_sparse/ # global presence index, sparse-packed (query benchmark)
|
||||
├── global_index_count/ # merged global count index
|
||||
├── specific_index_presence/ # species-specific presence indexes
|
||||
├── specific_index_count/ # species-specific count indexes
|
||||
├── query_dense/ # query output against global_index_presence
|
||||
├── query_sparse/ # query output against global_index_presence_sparse
|
||||
└── stats/ # all benchmark statistics
|
||||
├── indexing_presence/
|
||||
├── indexing_count/
|
||||
@@ -128,5 +172,9 @@ benchmark/
|
||||
├── specific_kmer_presence/
|
||||
├── specific_kmer_count/
|
||||
├── verify_merge_presence/
|
||||
└── verify_merge_count/
|
||||
├── verify_merge_count/
|
||||
├── pack_sparse/
|
||||
├── query_dense/
|
||||
├── query_sparse/
|
||||
└── verify_query/
|
||||
```
|
||||
|
||||
@@ -24,6 +24,12 @@ case "${TYPE}" in
|
||||
specific_kmer_presence|specific_kmer_count)
|
||||
HEADER="run,species,rebuild_wall_s,rebuild_rss_b,pack_wall_s,pack_rss_b,filter_total_wall_s,filter_total_rss_b,select_wall_s,select_rss_b,select_total_wall_s,select_total_rss_b"
|
||||
;;
|
||||
query_dense|query_sparse)
|
||||
HEADER="run,species,strain,query_wall_s,query_rss_b,total_wall_s,total_rss_b"
|
||||
;;
|
||||
verify_query)
|
||||
HEADER="run,species,strain,n_reads,n_common,missing_in_dense,missing_in_sparse,mismatched,mismatch_pct"
|
||||
;;
|
||||
*)
|
||||
echo "ERROR: unknown stats type '${TYPE}'" >&2
|
||||
exit 1
|
||||
|
||||
@@ -197,3 +197,17 @@ specific_index_count/Yersinia_ruckeri/index.done stats/specific_kmer_count/Yersi
|
||||
# Candidozyma_auris
|
||||
specific_index_presence/Candidozyma_auris/index.done stats/specific_kmer_presence/Candidozyma_auris.stats: global_index_presence/index.done
|
||||
specific_index_count/Candidozyma_auris/index.done stats/specific_kmer_count/Candidozyma_auris.stats: global_index_count/index.done
|
||||
|
||||
QUERY_SPECIMENS := Escherichia_coli--K-12_MG1655 Saccharolobus_islandicus--M.16.4
|
||||
|
||||
# query: Escherichia_coli--K-12_MG1655
|
||||
query_data/Escherichia_coli/K-12_MG1655/reads_R1.fastq.gz: genomes/GCF_000005845.2_ASM584v2_genomic.fna.gz
|
||||
query_dense/Escherichia_coli--K-12_MG1655.fasta.gz stats/query_dense/Escherichia_coli--K-12_MG1655.stats: query_data/Escherichia_coli/K-12_MG1655/reads_R1.fastq.gz global_index_presence/index.done
|
||||
query_sparse/Escherichia_coli--K-12_MG1655.fasta.gz stats/query_sparse/Escherichia_coli--K-12_MG1655.stats: query_data/Escherichia_coli/K-12_MG1655/reads_R1.fastq.gz global_index_presence_sparse/index.done
|
||||
stats/verify_query/Escherichia_coli--K-12_MG1655.stats: query_dense/Escherichia_coli--K-12_MG1655.fasta.gz query_sparse/Escherichia_coli--K-12_MG1655.fasta.gz
|
||||
|
||||
# query: Saccharolobus_islandicus--M.16.4
|
||||
query_data/Saccharolobus_islandicus/M.16.4/reads_R1.fastq.gz: genomes/GCF_000022445.1_ASM2244v1_genomic.fna.gz
|
||||
query_dense/Saccharolobus_islandicus--M.16.4.fasta.gz stats/query_dense/Saccharolobus_islandicus--M.16.4.stats: query_data/Saccharolobus_islandicus/M.16.4/reads_R1.fastq.gz global_index_presence/index.done
|
||||
query_sparse/Saccharolobus_islandicus--M.16.4.fasta.gz stats/query_sparse/Saccharolobus_islandicus--M.16.4.stats: query_data/Saccharolobus_islandicus/M.16.4/reads_R1.fastq.gz global_index_presence_sparse/index.done
|
||||
stats/verify_query/Saccharolobus_islandicus--M.16.4.stats: query_dense/Saccharolobus_islandicus--M.16.4.fasta.gz query_sparse/Saccharolobus_islandicus--M.16.4.fasta.gz
|
||||
|
||||
@@ -13,6 +13,11 @@ STOP_WORDS = {'complete', 'chromosome', 'whole', 'sequence', 'genome',
|
||||
'endosymbiont', 'of'}
|
||||
STOP_PREFIXES = ('scaffold', 'contig', 'plasmid')
|
||||
|
||||
# Specimens used as read sources for the query benchmark (see
|
||||
# DevDocMD/implementation/benchmark_query_testing.md): one common bacterium,
|
||||
# one distant lineage (the only archaeon in SPECIES).
|
||||
QUERY_SPECIMENS = ['Escherichia_coli--K-12_MG1655', 'Saccharolobus_islandicus--M.16.4']
|
||||
|
||||
|
||||
def is_stop(tok):
|
||||
t = tok.lower()
|
||||
@@ -113,6 +118,27 @@ def main():
|
||||
print(f'{sp_done} {sp_stats}: global_index_presence/index.done')
|
||||
print(f'{sc_done} {sc_stats}: global_index_count/index.done')
|
||||
|
||||
print()
|
||||
print('QUERY_SPECIMENS :=', ' '.join(QUERY_SPECIMENS))
|
||||
|
||||
by_specimen = {e[0]: e for e in entries}
|
||||
for specimen in QUERY_SPECIMENS:
|
||||
_, species, sim_dir, genome = by_specimen[specimen]
|
||||
query_dir = sim_dir.replace('simulated_data/', 'query_data/', 1)
|
||||
reads = f'{query_dir}/reads_R1.fastq.gz'
|
||||
dense_out = f'query_dense/{specimen}.fasta.gz'
|
||||
dense_stat = f'stats/query_dense/{specimen}.stats'
|
||||
sparse_out = f'query_sparse/{specimen}.fasta.gz'
|
||||
sparse_stat = f'stats/query_sparse/{specimen}.stats'
|
||||
vq_stat = f'stats/verify_query/{specimen}.stats'
|
||||
|
||||
print()
|
||||
print(f'# query: {specimen}')
|
||||
print(f'{reads}: {genome}')
|
||||
print(f'{dense_out} {dense_stat}: {reads} global_index_presence/index.done')
|
||||
print(f'{sparse_out} {sparse_stat}: {reads} global_index_presence_sparse/index.done')
|
||||
print(f'{vq_stat}: {dense_out} {sparse_out}')
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
|
||||
Executable
+86
@@ -0,0 +1,86 @@
|
||||
#!/usr/bin/env bash
|
||||
# Builds global_index_presence_sparse/ from global_index_presence/ by
|
||||
# copying the index (column files are kept in place after merge's dense
|
||||
# pack — see obikindex::KmerIndex::pack_matrices) and repacking in place
|
||||
# with --sparse.
|
||||
# Outputs:
|
||||
# global_index_presence_sparse/index.done (copied from source)
|
||||
# stats/pack_sparse/current.stats (one CSV data row, no header)
|
||||
set -euo pipefail
|
||||
|
||||
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
|
||||
BINARY="${SCRIPT_DIR}/../src/target/release/obikmer"
|
||||
SOURCE="${SCRIPT_DIR}/global_index_presence"
|
||||
OUTPUT="${SCRIPT_DIR}/global_index_presence_sparse"
|
||||
STATS_DIR="${SCRIPT_DIR}/stats/pack_sparse"
|
||||
STATS_FILE="${STATS_DIR}/current.stats"
|
||||
|
||||
mkdir -p "${STATS_DIR}"
|
||||
|
||||
echo "[pack_sparse] ${SOURCE} → ${OUTPUT}"
|
||||
|
||||
rm -rf "${OUTPUT}"
|
||||
cp -r "${SOURCE}" "${OUTPUT}"
|
||||
|
||||
STDERR_LOG=$(mktemp)
|
||||
trap 'rm -f "${STDERR_LOG}"' EXIT
|
||||
|
||||
"${BINARY}" pack --sparse "${OUTPUT}" 2>"${STDERR_LOG}"
|
||||
|
||||
cat "${STDERR_LOG}" >&2
|
||||
|
||||
python3 - "${STDERR_LOG}" <<'PYEOF' >"${STATS_FILE}"
|
||||
import sys, re
|
||||
|
||||
logfile = sys.argv[1]
|
||||
|
||||
def strip_ansi(s):
|
||||
return re.sub(r'\x1b\[[\x30-\x3f]*[\x20-\x2f]*[\x40-\x7e]', '', s)
|
||||
|
||||
def parse_wall(s):
|
||||
s = s.strip()
|
||||
if s.endswith('ms'): return float(s[:-2]) / 1000.0
|
||||
if s.endswith('s'): return float(s[:-1])
|
||||
return 0.0
|
||||
|
||||
def parse_rss(s):
|
||||
m = re.match(r'([\d.]+)\s*(GB|MB|KB|B)', s.strip())
|
||||
if not m: return 0
|
||||
return int(float(m.group(1)) * {'GB': 1<<30, 'MB': 1<<20, 'KB': 1024, 'B': 1}[m.group(2)])
|
||||
|
||||
def is_sep(s):
|
||||
return bool(s) and not re.search(r'[A-Za-z0-9]', s)
|
||||
|
||||
stats = {}
|
||||
state = 'scan'
|
||||
with open(logfile, errors='replace') as fh:
|
||||
for raw in fh:
|
||||
line = strip_ansi(raw.rstrip('\n'))
|
||||
s = line.strip()
|
||||
if state == 'scan':
|
||||
if re.search(r'\bstage\b.*\bwall\b', line):
|
||||
state = 'in_header'
|
||||
elif state == 'in_header':
|
||||
if is_sep(s): state = 'rows'
|
||||
elif state == 'rows':
|
||||
if is_sep(s): state = 'total'
|
||||
elif s:
|
||||
parts = re.split(r' +', s)
|
||||
if len(parts) >= 4:
|
||||
stats[parts[0]] = (parse_wall(parts[1]), parse_rss(parts[3]))
|
||||
elif state == 'total':
|
||||
if s:
|
||||
parts = re.split(r' +', s)
|
||||
if len(parts) >= 3:
|
||||
stats['TOTAL'] = (parse_wall(parts[1]),
|
||||
parse_rss(parts[3]) if len(parts) > 3 else 0)
|
||||
break
|
||||
|
||||
w, r = stats.get('pack', ('', ''))
|
||||
tw, tr = stats.get('TOTAL', ('', ''))
|
||||
row = [f'{w:.3f}' if isinstance(w, float) else '', str(r),
|
||||
f'{tw:.3f}' if isinstance(tw, float) else '', str(tr)]
|
||||
print(','.join(row))
|
||||
PYEOF
|
||||
|
||||
echo "Done → ${OUTPUT}"
|
||||
Executable
+105
@@ -0,0 +1,105 @@
|
||||
#!/usr/bin/env bash
|
||||
# Usage: query_one.sh MODE SPECIMEN
|
||||
# MODE = dense | sparse
|
||||
# SPECIMEN = "species--strain" (Make pattern stem), reads from query_data/
|
||||
# Outputs:
|
||||
# query_MODE/SPECIMEN.fasta.gz (obikmer query output, --count-missing)
|
||||
# stats/query_MODE/SPECIMEN.stats (one CSV data row, no header)
|
||||
set -euo pipefail
|
||||
|
||||
MODE="$1"
|
||||
SPECIMEN="$2"
|
||||
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
|
||||
BINARY="${SCRIPT_DIR}/../src/target/release/obikmer"
|
||||
|
||||
case "${MODE}" in
|
||||
dense) INDEX="${SCRIPT_DIR}/global_index_presence" ;;
|
||||
sparse) INDEX="${SCRIPT_DIR}/global_index_presence_sparse" ;;
|
||||
*) echo "ERROR: unknown mode '${MODE}' (expected dense|sparse)" >&2; exit 1 ;;
|
||||
esac
|
||||
|
||||
species="${SPECIMEN%%--*}"
|
||||
strain="${SPECIMEN#*--}"
|
||||
|
||||
READS_DIR="${SCRIPT_DIR}/query_data/${species}/${strain}"
|
||||
OUT_DIR="${SCRIPT_DIR}/query_${MODE}"
|
||||
STATS_DIR="${SCRIPT_DIR}/stats/query_${MODE}"
|
||||
OUT_FILE="${OUT_DIR}/${SPECIMEN}.fasta.gz"
|
||||
STATS_FILE="${STATS_DIR}/${SPECIMEN}.stats"
|
||||
|
||||
mkdir -p "${OUT_DIR}" "${STATS_DIR}"
|
||||
|
||||
r1="${READS_DIR}/reads_R1.fastq.gz"
|
||||
r2="${READS_DIR}/reads_R2.fastq.gz"
|
||||
if [[ ! -f "${r1}" || ! -f "${r2}" ]]; then
|
||||
echo "ERROR: reads not found in ${READS_DIR}" >&2
|
||||
exit 1
|
||||
fi
|
||||
|
||||
echo "[${SPECIMEN}] query (${MODE}) → ${OUT_FILE}"
|
||||
|
||||
STDERR_LOG=$(mktemp)
|
||||
trap 'rm -f "${STDERR_LOG}"' EXIT
|
||||
|
||||
"${BINARY}" query \
|
||||
--count-missing \
|
||||
"${INDEX}" "${r1}" "${r2}" \
|
||||
2>"${STDERR_LOG}" \
|
||||
| gzip >"${OUT_FILE}"
|
||||
|
||||
cat "${STDERR_LOG}" >&2
|
||||
|
||||
python3 - "${species}" "${strain}" "${STDERR_LOG}" <<'PYEOF' >"${STATS_FILE}"
|
||||
import sys, re
|
||||
|
||||
species, strain, logfile = sys.argv[1], sys.argv[2], sys.argv[3]
|
||||
|
||||
def strip_ansi(s):
|
||||
return re.sub(r'\x1b\[[\x30-\x3f]*[\x20-\x2f]*[\x40-\x7e]', '', s)
|
||||
|
||||
def parse_wall(s):
|
||||
s = s.strip()
|
||||
if s.endswith('ms'): return float(s[:-2]) / 1000.0
|
||||
if s.endswith('s'): return float(s[:-1])
|
||||
return 0.0
|
||||
|
||||
def parse_rss(s):
|
||||
m = re.match(r'([\d.]+)\s*(GB|MB|KB|B)', s.strip())
|
||||
if not m: return 0
|
||||
return int(float(m.group(1)) * {'GB': 1<<30, 'MB': 1<<20, 'KB': 1024, 'B': 1}[m.group(2)])
|
||||
|
||||
def is_sep(s):
|
||||
return bool(s) and not re.search(r'[A-Za-z0-9]', s)
|
||||
|
||||
stats = {}
|
||||
state = 'scan'
|
||||
with open(logfile, errors='replace') as fh:
|
||||
for raw in fh:
|
||||
line = strip_ansi(raw.rstrip('\n'))
|
||||
s = line.strip()
|
||||
if state == 'scan':
|
||||
if re.search(r'\bstage\b.*\bwall\b', line):
|
||||
state = 'in_header'
|
||||
elif state == 'in_header':
|
||||
if is_sep(s): state = 'rows'
|
||||
elif state == 'rows':
|
||||
if is_sep(s): state = 'total'
|
||||
elif s:
|
||||
parts = re.split(r' +', s)
|
||||
if len(parts) >= 4:
|
||||
stats[parts[0]] = (parse_wall(parts[1]), parse_rss(parts[3]))
|
||||
elif state == 'total':
|
||||
if s:
|
||||
parts = re.split(r' +', s)
|
||||
if len(parts) >= 3:
|
||||
stats['TOTAL'] = (parse_wall(parts[1]),
|
||||
parse_rss(parts[3]) if len(parts) > 3 else 0)
|
||||
break
|
||||
|
||||
qw, qr = stats.get('query', ('', ''))
|
||||
tw, tr = stats.get('TOTAL', ('', ''))
|
||||
row = [species, strain,
|
||||
f'{qw:.3f}' if isinstance(qw, float) else '', str(qr),
|
||||
f'{tw:.3f}' if isinstance(tw, float) else '', str(tr)]
|
||||
print(','.join(row))
|
||||
PYEOF
|
||||
Executable
+34
@@ -0,0 +1,34 @@
|
||||
#!/usr/bin/env bash
|
||||
# Usage: simulate_query_one.sh genome.fna.gz output_dir
|
||||
# Simulates a fixed-size paired-end HiSeq read set for the query benchmark.
|
||||
# Unlike simulate_one.sh (coverage-proportional, used to build the indexed
|
||||
# specimens), this always draws N_READS pairs regardless of genome size —
|
||||
# query benchmark numbers (wall/RSS) must stay comparable across genomes of
|
||||
# very different sizes. Independent iss run (unseeded), so error draw
|
||||
# differs from any simulated_data/ reads for the same genome.
|
||||
set -euo pipefail
|
||||
|
||||
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
|
||||
ISS="${SCRIPT_DIR}/../.venv/bin/iss"
|
||||
N_READS=100000
|
||||
CPUS="${CPUS:-$(sysctl -n hw.logicalcpu 2>/dev/null || nproc 2>/dev/null || echo 2)}"
|
||||
|
||||
genome_file="$1"
|
||||
out_dir="$2"
|
||||
|
||||
mkdir -p "${out_dir}"
|
||||
|
||||
tmp_fasta=$(mktemp "${TMPDIR:-/tmp}/obikmer_XXXXXX.fna")
|
||||
trap 'rm -f "${tmp_fasta}"' EXIT
|
||||
|
||||
gzip -dc "${genome_file}" > "${tmp_fasta}"
|
||||
|
||||
echo "[${out_dir}] ${N_READS} read pairs (query benchmark, fixed size)"
|
||||
|
||||
"${ISS}" generate \
|
||||
--genomes "${tmp_fasta}" \
|
||||
--model HiSeq \
|
||||
--n_reads "${N_READS}" \
|
||||
--cpus "${CPUS}" \
|
||||
--compress \
|
||||
--output "${out_dir}/reads"
|
||||
Executable
+68
@@ -0,0 +1,68 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Compare dense vs sparse obikmer query output for one specimen.
|
||||
|
||||
Both files are `obikmer query --count-missing` output: OBITools4-style
|
||||
FASTA, one record per read — `>id {"kmer_count":N,"kmer_missing":M,
|
||||
"kmer_strict_matches":{"label":count,...}}`. Packing format (dense vs
|
||||
sparse presence matrix) must not change query results — only I/O access
|
||||
pattern differs. Matched by read id rather than by stream position: the
|
||||
query pipeline processes input in chunks across worker threads and does
|
||||
not guarantee output order matches input order.
|
||||
|
||||
Output to stdout: one CSV row
|
||||
species, strain, n_reads, n_common, missing_in_dense, missing_in_sparse, mismatched, mismatch_pct
|
||||
"""
|
||||
import argparse
|
||||
import gzip
|
||||
import json
|
||||
import sys
|
||||
|
||||
|
||||
def load(path):
|
||||
records = {}
|
||||
opener = gzip.open if path.endswith('.gz') else open
|
||||
with opener(path, 'rt') as fh:
|
||||
for line in fh:
|
||||
if not line.startswith('>'):
|
||||
continue
|
||||
header = line[1:].rstrip('\n')
|
||||
read_id, _, json_part = header.partition(' ')
|
||||
records[read_id] = json.loads(json_part) if json_part else {}
|
||||
return records
|
||||
|
||||
|
||||
def annotations_equal(a, b):
|
||||
return (
|
||||
a.get('kmer_count') == b.get('kmer_count')
|
||||
and a.get('kmer_missing') == b.get('kmer_missing')
|
||||
and a.get('kmer_strict_matches', {}) == b.get('kmer_strict_matches', {})
|
||||
)
|
||||
|
||||
|
||||
def main():
|
||||
ap = argparse.ArgumentParser()
|
||||
ap.add_argument('--species', required=True)
|
||||
ap.add_argument('--strain', required=True)
|
||||
ap.add_argument('dense_fasta')
|
||||
ap.add_argument('sparse_fasta')
|
||||
args = ap.parse_args()
|
||||
|
||||
dense = load(args.dense_fasta)
|
||||
sparse = load(args.sparse_fasta)
|
||||
|
||||
dense_ids, sparse_ids = set(dense), set(sparse)
|
||||
common = dense_ids & sparse_ids
|
||||
missing_in_dense = len(sparse_ids - dense_ids)
|
||||
missing_in_sparse = len(dense_ids - sparse_ids)
|
||||
|
||||
mismatched = sum(1 for rid in common if not annotations_equal(dense[rid], sparse[rid]))
|
||||
|
||||
n_reads = len(dense_ids | sparse_ids)
|
||||
mismatch_pct = 100.0 * (mismatched + missing_in_dense + missing_in_sparse) / n_reads if n_reads else 0.0
|
||||
|
||||
print(f'{args.species},{args.strain},{n_reads},{len(common)},'
|
||||
f'{missing_in_dense},{missing_in_sparse},{mismatched},{mismatch_pct:.6f}')
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
Executable
+28
@@ -0,0 +1,28 @@
|
||||
#!/usr/bin/env bash
|
||||
# Usage: verify_query_one.sh SPECIMEN
|
||||
# SPECIMEN = "species--strain" (Make pattern stem)
|
||||
# Output: stats/verify_query/SPECIMEN.stats (one CSV data row, no header)
|
||||
set -euo pipefail
|
||||
|
||||
SPECIMEN="$1"
|
||||
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
|
||||
PYTHON="${SCRIPT_DIR}/../.venv/bin/python3"
|
||||
VERIFY_PY="${SCRIPT_DIR}/verify_query.py"
|
||||
|
||||
species="${SPECIMEN%%--*}"
|
||||
strain="${SPECIMEN#*--}"
|
||||
|
||||
DENSE="${SCRIPT_DIR}/query_dense/${SPECIMEN}.fasta.gz"
|
||||
SPARSE="${SCRIPT_DIR}/query_sparse/${SPECIMEN}.fasta.gz"
|
||||
STATS_DIR="${SCRIPT_DIR}/stats/verify_query"
|
||||
STATS_FILE="${STATS_DIR}/${SPECIMEN}.stats"
|
||||
|
||||
mkdir -p "${STATS_DIR}"
|
||||
|
||||
echo "[${SPECIMEN}] verifying query (dense vs sparse)"
|
||||
|
||||
"${PYTHON}" "${VERIFY_PY}" \
|
||||
--species "${species}" \
|
||||
--strain "${strain}" \
|
||||
"${DENSE}" "${SPARSE}" \
|
||||
>"${STATS_FILE}"
|
||||
+1
-1
@@ -55,7 +55,7 @@ nav:
|
||||
- Kmer filtering: implementation/filtering.md
|
||||
- Select command: implementation/select.md
|
||||
- obitaxonomy crate: implementation/obitaxonomy.md
|
||||
- "Benchmark: query-path testing (discussion)": implementation/benchmark_query_testing.md
|
||||
- "Benchmark: query-path testing": implementation/benchmark_query_testing.md
|
||||
- Architecture:
|
||||
- Sequences: architecture/sequences/invariant.md
|
||||
- Kmer index: architecture/index_architecture.md
|
||||
|
||||
Reference in New Issue
Block a user