BrucellaID 2.2.0

The extended reference set changes, so some samples change their call. It is now 429 genomes instead of 432, and sixteen of them carry a different label. Distances do not change: a sample's distance to a reference does not depend on what else the sketch holds, and the 429 remaining genomes are in the rebuilt sketch with exactly the hashes they had.

A sample can now be reported under a species name that does not begin with Brucella. Ochrobactrum teleogrylli and O. chromiisoli were described after Ochrobactrum was moved into Brucella, nobody has published a Brucella combination for either, and NCBI files both outside the genus. Four genomes of the set are theirs, at 98.1–99.2% ANI to the type strains, and are labelled with the names the species have. Anything that parses the Species column on the assumption that it starts with "Brucella" should be checked.

Eleven unnamed genomes take the name of the Brucella species they belong to. They were Brucella sp. and are 96.2–99.3% ANI from a type strain, measured in both directions: one B. tritici, three B. anthropi, two B. intermedia, one B. haematophila and four B. pseudogrignonensis. A sample matching one was reported as Brucella sp. and is now reported under its species.

GCA_900095155.1 is now Brucella sp., not B. inopinata. It is 97.87% ANI from the B. inopinata type strain BO1, further than three references the set already called Brucella sp.

Three partial genomes are removed: GCA_040391195.1, GCA_032085995.1 and GCA_032104425.1, each 2.1 to 3.7 Mb of a species of about 5 Mb that the set holds complete genomes of. None was the best match of any identified run among 4,024 public read sets. A custom --references table written against the 432 will name three genomes the sketch no longer holds, which BrucellaID reports at start-up.

What it changes in practice

Population	Calls that change
The 4,024 usable public Brucella read sets in ENA	four: SRR12605656, SRR25828924 and SRR25828932 go from Brucella sp. to O. teleogrylli, and SRR12224742 from B. inopinata to Brucella sp.
1,543 held-out simulated samples of a named species	none

The 16-genome vSNP set (--reference-set vsnp) is untouched.

Arbitration also reaches fewer references. Genome against genome, a second species falls inside the close-match margin of 33 of the 429 references across 10 species pairs, where it was 38 across 14: four of those pairs were a named species against one of its own unnamed genomes.

The evidence for every relabel and removal is in curation/unnamed_lineages.md, and each decision is a row of scripts/extended_references_curation.tsv.

Also in this release

• The at-scale validation is re-measured under the shipped k-mer copy threshold. All 4,024 public read sets were swept again against both reference sets. The headline figures are unchanged: the extended set matches the submitter on 98.7% and the vSNP set on 96.2%.
• --loose-distance 0.001 is confirmed on all 4,001 identified runs, where the 95th percentile of agreeing calls is 0.00068.
• One claim is withdrawn. Agreement with the submitter is not flat in distance: it falls from 99.5% on close calls to 81.2% between 0.005 and 0.02.
• A warning for ENA's three-file runs. find_samples() now warns when a _3 fastq sits beside a pair it has formed: in that layout _1 holds the unpaired reads and _2/_3 are the mates.
• A stale sentence on the validation page is corrected: of the 95 unnamed simulated samples, 37 are reported as Brucella sp. and 56 under a name.

The full entry is in the changelog.

BrucellaID 2.1.1

The quality assurance report no longer shows the biovars recorded for a sample's sequence type. The
Biovars of its isolates row is gone from each sample's evidence block, in both the HTML and the PDF.

The reason is the shape of the data rather than the idea. Of the 118 sequence types in the shipped table,
42 have a single biovar on record and 59 have none, so on about a third of samples the row rendered one
bare value — 2:60 — beside that sample's own species and sequence type. That reads as a biovar call,
which is how the reference-genome biovar column came to be misread before 1.1.0 removed it. The 17 sequence
types that record several biovars did convey the ambiguity well, but they were the minority.

Nothing measured is lost. The counts stay in the packaged mlst_st_species.tsv, in a custom
--st-species table, and in brucellaid.json. The reference genomes' own biovars stay in
brucellaid_mash.tsv, where they are facts about the references rather than calls about the sample. The
report's methods section still sets out why no biovar can be measured — 7 of the 432 references carry one,
the two B. suis biovar 2 genomes sit further apart than the biovar 3 genome is from one of them, and the
best any sequence-type-to-biovar rule manages over the isolate database is 66% — and now says plainly that
this is why it shows none at all.

Also in this release:

• The wiki's biovar explanation is cut to one short treatment under
  How it works, where it had grown
  to seven places and about forty lines of argument for something the tool does not do.
• 393 tests, 100% statement and branch coverage. One is new, and asserts the row is absent from both
  reports while the fixture still carries biovars, so it cannot come back unnoticed.

No identification changes: the species call, the sequence type, the distances and the shared-hash counts
are all as in 2.1.0.

BrucellaID 2.1.0

The k-mer copy threshold was re-derived, and reported distances change because of it. Before MASH
sketches a read set it drops every k-mer seen fewer than m times, to throw away the ones that exist only
because of a sequencing error. m was one more copy per 20x of coverage, capped at 20; it is now one more
copy per 5x, capped at 48. The floor of 2 is unchanged.

No species call and no sequence type changes: the 2,022-sample simulated validation returns exactly the
same 1,538 of 1,543 held-out named samples, and the repository's own suite still passes 22 of 22. What
changes is MashDistance and SharedHashes, on real libraries especially:

library	depth	before	after
ERR2136546	47x	0.00652, 773/1000	0.000119, 995/1000
ERR3454415	70x	0.00068, 972/1000	0.000144, 994/1000
SRR19520329	53x	0.00051, 979/1000	0.000119, 995/1000
SRR8550484	69x	0.00645, 775/1000	0.000264, 989/1000

If you compare distances across versions, they are not comparable. A distance from 2.0.0 or earlier is
larger than the same measurement today, by a factor that depends on the depth of the library. Pin any fixed
value with --min-kmer-copies to reproduce an old report: 2 gives a 1.1.0 report exactly.

What it was fitted on. 81 threshold sweeps over 63 real libraries from 5x to 1903x. The new rule lands
inside the measured plateau for 80 of the 81 curves; the old one for 26. Simulated reads cannot settle it,
and an earlier derivation that rested on them got it wrong.

Also in this release

• --loose-distance is recalibrated, 0.005 to 0.001. Same rule as before — the 95th percentile of the
  distance of calls that agree with the ENA submitter — but that percentile was measured on distances the
  old copy threshold had inflated about sevenfold. Re-measured on 207 public read sets it is 0.00066. The
  note stays about as rare, flagging 4.0% of identified runs against 4.9% before.
• A claim behind the old calibration was circular, and is withdrawn. It offered as corroboration that
  every agreeing call past the threshold shared fewer than 900 of its 1000 hashes. MASH derives the
  distance from the shared-hash count, so that is the threshold restated in another unit. What does
  support leaving the species call alone is that agreement with the submitter is flat in distance.
• Three head-to-head figures of the public sweep are corrected, all stale the same way and all moving
  in the extended set's favour: it wins 94 of the 146 disagreements rather than 91, 89 of them are runs
  the vSNP set cannot name rather than 72, and 21 of those are B. intermedia rather than 15.
• The loose-match note's advice is reversed where it suggested raising --min-kmer-copies; on a deep
  mixed library the threshold now has to come down.
• Evidence is shipped for every derived constant: scripts/kmer_threshold_sweep.tsv,
  kmer_threshold_simulated.tsv, loose_distance_calibration.tsv,
  public_sweep_disagreements.tsv, with recalibrate_loose_distance.py and score_public_sweep.py
  deriving the published figures from them.
• A fix in validation/check_results.py, which after the rule change could fail a sample whose coverage
  it had just accepted.
• 392 tests, 100% statement and branch coverage.

Full detail: Changelog ·
The k-mer copy threshold

BrucellaID 2.0.0

The extended reference set is now the default, and that is the breaking change. A run with no --reference-set identifies against 432 curated genomes covering all 30 named species of the Brucella clade — the former Ochrobactrum included — instead of vSNP's 16 genomes and 8 species. Over every usable public Brucella read set in ENA, 4,042 runs, it matches the submitter on 98.7% against the 16-genome set's 96.2%; the two disagree on 146 runs, where the extended set is right on 91 of them and the smaller set on none. Across the 2,018 B. melitensis and 1,118 B. abortus runs that carry most real traffic the two agree, so the change adds coverage without disturbing the calls the smaller set already made well. The numbers are in the wiki.

What this costs a 1.x pipeline. No table changed shape — not one column was added or removed anywhere — but three things a parser might read did change:

• ClosestReference is a GCA file name now (GCA_000742195.1_ASM74219v1_genomic.fna.gz), and MashDistance and SharedHashes are measured against a different, usually closer, genome.
• The code column of --xlsx and species_code in the JSON hold a species epithet (abortus, abortus5) rather than a vSNP code. The break is partial, which is worse than total: 7 of vSNP's 16 codes — canis, microti, neotomae, ovis, suis1, suis2, suis3 — spell the same thing in both tables and keep matching, while ab1, ab3, mel1, mel1b, mel2, mel3, ceti1, ceti2 and suis4 silently stop. A code-to-species lookup will half work, per species. --xlsx now warns when it is writing epithets.
• Notes is populated on nearly every sample, where in 1.x it was usually empty, because the dense set almost always has another genome of the same species inside the close-match margin. On the four tutorial read sets it goes from 1 of 4 to 4 of 4. A rule of "a non-empty Notes means look at this one" now fires on everything.

brucellaid_mash.tsv holds every distance, so it grows with the set: 5 kB to 170 kB for four samples, and brucellaid.json 24 kB to 132 kB.

--reference-set vsnp gives back 1.x's reference data — the same 16 genomes, the same sketch, the same threshold, and therefore the same species, closest reference, distance, shared hashes and sequence type on every sample we have checked. It is not a byte-for-byte 1.x run: the 2.0 features that do not depend on the reference set still apply, so a library past --loose-distance gains a note, inserted first, and Loose match beyond appears in the run parameters. --reference-set default, which only ever existed in development builds, is rejected with exit 2 rather than quietly accepted.

The default moved at the library level too, not only on the command line: reference_sketch() and load_references() resolve to the extended files when given none, so a script that calls analyse_samples() directly measures against the 432 genomes without knowing the flag exists.

Two other things are new. A closest reference that is not close at all now says so: past --loose-distance, 0.005 by default, the call carries a note with the distance and the shared hashes, without changing the species. The threshold comes from the ENA sweep, where nineteen in twenty calls that agreed with the submitter sat within 0.005. And a nearly-as-close reference of the same species is counted rather than listed, so the note reads "4 other Brucella melitensis genome(s)" instead of four file names that change nothing.

One of the 30 species took some finding: nothing in INSDC is labelled B. papionis, but its genome is there under "Brucella sp.", and how we identified it is in the repository, along with every mislabelled public record the validation turned up, all of them reported to the archives.

Fixed while reviewing this release, none of it a consequence of the new default: the legacy --xlsx report wrote a value beginning with = as a live Excel formula, so a sample name could display something other than the sample; a named pipe or a link to a device in an input folder could hang a run or end it before any report was written, losing every other sample; and a species resting on the sequence type alone, with no reference behind it, now counts as a warning rather than "ok" in the quality assurance report.

One caveat worth knowing, which predates 2.0 and is unchanged by it: the derived k-mer copy threshold is still low for some libraries, and on two of the four tutorial read sets it leaves MashDistance inflated roughly 25-fold — --min-kmer-copies 12 takes SRR8550484 from 775/1000 hashes at 0.0065 to 990/1000 at 0.00024. The species and the sequence type are right either way, and --reference-set vsnp gives the same figures, so read a distance from a default run as a floor on the match rather than a measurement of it. The evidence is on the Validation page.

See an example report, and the full list is in the changelog.

BrucellaID 1.2.1

A patch release: what a full review of the code and an audit of the documentation caught after 1.2.0.

Concatenated gzip files are sized honestly. The coverage estimate trusted the gzip trailer, which names only the last member of a file — so lane files joined with cat read as a fraction of their size, and a bgzf file (what bgzip and every htslib tool write) as almost nothing. A deep library could estimate under 1x, keep the k-mer copy threshold at its floor of 2, and let the error k-mers back into the sketch: the very failure 1.2.0 exists to prevent. The head of the file is now decompressed member by member, several members mean the trailer is ignored, and a trailer that contradicts the measured compression ratio is ignored too. This is the fix to have if you feed BrucellaID merged or bgzipped fastq files.

A run is harder to lose. The TSV, JSON and MultiQC reports are written in UTF-8 whatever the locale, so a sample name with an accent no longer costs a finished run its reports under LC_ALL=C; and a gzip file with a corrupt stream inside fails that one sample with the real error instead of crashing the whole --parallel run.

The reports read better. Tables are only as wide as their content instead of being stretched across the page, every column is left-aligned, and a rule opens and closes each table — in the HTML and the PDF alike. A coverage below 10x is written with one decimal ("7.5x", "under 0.1x"), never "roughly 0x". The coverage log line names the sample it describes, so parallel samples stay apart. See the example report.

The documentation agrees with the code again, down to the methods text inside the generated report: the last biovar ghosts are gone (BrucellaID has reported none since 1.1.0), the k-mer filter is described as the depth-derived threshold everywhere, and Output-files documents min_kmer_copies and estimated_coverage.

The full list is in the changelog.

BrucellaID 1.2.0

The k-mer copy threshold now rises with the sequencing depth. mash dist -m 2 keeps every k-mer a read set holds twice — and an error k-mer turns up about depth × error rate times, so past roughly 50x the errors reach two copies on their own and crowd the true k-mers out of the bottom-1000 sketch. The distance then grows with the depth instead of with the genetic distance, and above about 240x the sample got no call at all. The threshold is now derived from each sample — one more required copy per 20x of coverage, within [2, 20] — with the coverage estimated from the read files in a tenth of a second, without decompressing them. Below about 50x nothing changes. --min-kmer-copies N pins any fixed value, including the old 2.

Validated at depth, and on real reads. The validation now simulates two genomes at 150x and 250x, where the derived threshold keeps 1000/1000 sketch hashes and the old fixed 2 kept 728–902, and runs the four public ENA libraries of the tutorial next to the simulated samples: each must give the species its submitter gave it, the sequence type it is known to carry, and the same reference as the simulated sample of its own reference genome. All 22 pass — and the fix shows on real data: SRR8550484 (~69x) moves from 422/1000 shared hashes at distance 0.0248 to 775/1000 at 0.0065. The numbers are in the wiki.

See an example report — it now says which threshold each sample got and the coverage it came from, as do brucellaid.json (min_kmer_copies, estimated_coverage) and the PDF.

Upgrading: no report or table changed shape. Results change only on read sets above about 50x, where they were the ones being distorted: distances shrink toward what an assembly would give, shared hashes rise, and a sample deep enough to have been ND can get its call back.

The full list is in the changelog.

BrucellaID 1.1.0

The sequence type now decides the species the MASH distance cannot. B. canis ATCC 23365 and B. suis bv. 4 strain 40 are 0.00024 apart in the sketch — the closest pair of the sixteen, closer than two B. suis references are to each other — so which one a canis read set lands on first is noise. The two species share no sequence type at all, so where MASH leaves both within reach, the ST chooses, from a new table of the species of every ST counted over 3,426 PubMLST isolates.

Every run writes a quality assurance report. brucellaid_report.html is one self-contained file — no script, no stylesheet, no font, nothing fetched from anywhere — holding the closest reference genomes, how the species was decided, every allele call with the margin it had, the input files with their MD5, the exact command, and the checksum of every reference data file the run measured against, with a box for a reviewer to sign. One sample to a printed page. --pdf writes the same report as a PDF, and brucellaid_mqc.json folds the run into MultiQC.

See an example report, or the same thing as a PDF.

The biovar column is gone, because BrucellaID cannot measure a biovar: the sketch has references for only B. abortus biovars 1 and 3, the two B. melitensis biovar 1 references are further apart than the biovar 2 and 3 references are from each other, and the best any ST → biovar rule could manage is 66% for B. abortus. The numbers are in the wiki.

Upgrading: brucellaid.tsv gained SpeciesSource and STSpecies and lost Biovar, so anything reading it by column number needs updating; by column name, only Biovar disappears. --no-st-species restores the 1.0.0 species call exactly.

The full list is in the changelog.

BrucellaID 1.0.0

BrucellaID becomes an installable, tested package. It answers the same two questions as the first version —
which Brucella species a read set is, and which sequence type — and answers the first one exactly as before.

Install

conda create -n brucellaid -c conda-forge -c bioconda brucellaid   # once the bioconda recipe is merged
# or
pip install brucellaid                                              # MASH and BBTools from bioconda

What changed

• The two OLC dependencies are gone: accessoryFunctions (OLCTools) and MLSTsippr (sipprverse). Neither
  installs from PyPI, so neither could be packaged. BrucellaID now needs nothing but the Python standard
  library, plus MASH and BBTools.
• The MLST database is in the package, so -d is gone; the input is -r1/-r2 or -i, instead of -s.
• The MLST alleles are called with Seal from BBTools, and the database is PubMLST's: 309 alleles and 122
  sequence types, where vSNP's copy had 263 and 81. That copy had drifted — it still held a gyrB_23 that
  PubMLST has deleted, and that sequence matched the B. abortus 9-941 genome better than the real gyrB_1,
  turning ST 1 into a new profile.
• The species identification is unchanged in substance: the same 16 vSNP reference genomes, the same sketch,
  the same threshold of 150 shared hashes out of 1000. New around it: the distance to every reference, the
  species/biovar/strain/accessions of each, and a note when a second reference was nearly as close.
• Four reports — one row per sample, one per locus, every MASH distance, and a JSON report — and --xlsx
  still writes stat_alignment_summary_<date>.xlsx for whatever reads it.
• 194 tests at 100% coverage, on Python 3.10, 3.12 and 3.14; a wiki; a tutorial on four public read sets.
• Validated on 16 genomes with a known answer: every species, biovar and closest reference right, all 14
  sequence types the one the genome carries, and the nearest neighbouring species (B. anthropi) reported
  as not identified at 4 shared hashes of 1000.

The 193 MB of per-reference vSNP script_dependents folders are no longer in the repository; they remain in
the history and in the v0.1 release.

Full changelog: https://github.com/duceppemo/BrucellaID/wiki/Changelog

BrucellaID 0.1

The clone-and-run version of BrucellaID: brucellaID.py and brucellaID_methods.py, run from a clone of the
repository, identifying the Brucella species of a read set with MASH against the 16 vSNP reference genomes and
typing it with the Whatmore MLST-9 scheme through MLSTsippr.

Tagged so that analyses and scripts built on this layout keep a fixed point to refer to: everything from
v1.0.0 onwards is an installed package with a
different command line.

Running it

git clone https://github.com/duceppemo/BrucellaID
cd BrucellaID
python brucellaID.py -s /path/to/fastq_folder -d /path/to/mlst_database -o /path/to/output

The MLST database -d points at is downloaded separately with get_mlst.py (from SRST2), which fetches the
Brucella scheme from PubMLST.

What it needs

• Python 3 with OLCTools (accessoryFunctions) and
  sipprverse (MLSTsippr)
• MASH and BBTools on the PATH

What it writes

stat_alignment_summary_<date>.xlsx in the output folder: one row per strain with the closest reference genome,
the species, and the MLST sequence type.

This release also contains the dependencies/ folder — the 16 reference genomes, the MASH sketch
vsnp_reference.msh, and the 193 MB of per-reference vSNP script_dependents trees — which later versions no
longer carry in the working tree.