TEsorter2

Fast, divergence-robust classification of transposable elements.

TEsorter2 is a reimplementation of TEsorter that keeps its classification semantics while introducing three major improvements:

1. Speed: HMMER's limited parallelism is replaced by pyHMMER with workload-aware load balancing, plus an optional facet pre-screen and a parallelized second-pass BLAST.
2. Sensitivity on degraded copies: an optional BATH engine performs frameshift-aware translated search directly against nucleotide sequence.
3. Reproducibility: clade assignment uses a score-weighted vote instead of a count-based vote whose ties were broken by internal data-structure ordering. It also adds multi-database reconciliation in a single run and a genome mode for both engines.

Choosing an engine

	Element mode (pre-extracted TEs)	Genome mode (assembly)
Intact / low-divergence sequence	pyHMMER (default) : fastest	pyHMMER or BATH
Degraded / frameshifted copies	BATH (--bath) : slower than pyHMMER, more sensitive	BATH (--bath) : both faster and more sensitive

---

Search modes

Default mode (pyHMMER)

Single-pass --nobias hmmsearch against all models via pyHMMER, in-process.

• Model-cost-aware parallel load balancing chooses between pyHMMER's queries and targets parallelization per model bin, reaching near-full CPU utilization (see Parallel strategy).
• Results are written to a SQLite database, so filtering and re-analysis do not require re-running the search.

Facet mode (--facet)

Pre-screens amino-acid databases with spliced sub-HMMs ("facets") to route each sequence only to the models likely to produce its best hit:

1. Facet screen: tiered sub-HMMs (96 → 64 → 48 → 32 aa) searched against all six translated frames.
2. Targeted verification: top facet hit per domain family verified with a full-model --nobias search.
3. Cross-family completion: verified frames searched for missing domain families.
4. Legacy fallback: frames with no facet signal get a full search. DNA databases always use the default search; DNA facets do not repay their overhead. Incompatible with --bath and --genome.

BATH mode (--bath)

Replaces pyHMMER with bathsearch --fs for amino-acid databases. BATH aligns protein pHMMs directly against raw nucleotide input, allowing the alignment to change frame at indels and read through stop codons under a frameshift penalty, so a domain interrupted by a frameshift is recovered as a single hit with its true extent. No six-frame translation is performed.

• HMMER3 databases are converted on first use (bathconvert) and cached as {db}.bath.hmm.
• Hits are normalized into the same internal schema as the HMMER path, so classification, reconciliation and BLAST pass-2 are unchanged.
• BATH tblout reports no per-residue posterior probability, so acc is set to 1.0 and the minimum-accuracy filter is a no-op for BATH hits.
• Minus-strand coordinates are normalized to ascending order, with strand encoded in the target suffix, matching the HMMER convention.
• Incompatible with --facet.

Genome mode (--genome)

Treats the input as whole-genome sequence rather than pre-extracted elements: detects TE protein domains throughout, classifies each domain individually, resolves overlapping features, and emits a domain-level GFF3 plus a summary table. It does not produce a per-element .cls.tsv and does not run BLAST pass-2, matching TEsorter's -genome behaviour.

• pyHMMER: six-frame translates each window, then maps amino-acid envelope coordinates back to nucleotide space. Window size is capped automatically to stay under pyHMMER's 100k-residue per-sequence limit; the retained overlap exceeds any TE domain, so no domains are lost.
• --bath: runs bathsearch --fs directly on nucleotide windows. The tblout already reports nucleotide coordinates and strand, so no translation or coordinate back-mapping is needed. BATH additionally streams long targets in ~0.25 Mb blocks overlapped by the maximum expected hit length, reconciling boundary duplicates internally. Window size: --win-size (default 1e6), --win-ovl (default 1e5). Requires at least one amino-acid database; DNA-only databases (sine) are skipped. Incompatible with --facet.

Multiple databases

TEsorter2 accepts several databases in one run (-d rexdb,gydb) and reconciles them in three stages:

1. Independent classification: each database classifies every element on its own, emitting a per-database {prefix}.{db}.cls.tsv in native TEsorter format.
2. Name harmonization: per-database calls are projected onto a unified taxonomy that collapses superfamily synonyms (e.g. Pao → Bel-Pao) and unifies clade names. Lineages with no established equivalent are kept distinct to avoid spurious agreement.
3. Scope-aware hierarchical reconciliation: a hierarchical vote at Order → Superfamily → Clade. At each level candidate labels are weighted by the summed normalized domain score per database, and only entries consistent with the winning label advance. The superfamily level is scope-aware: a database may only elect a superfamily it models with at least 2 clades. Every per-database call is retained in the SecondaryHits column as db:order/superfamily/clade=score, in descending order of evidence. Use --compat-tesorter-output for the original 7-column format.

Clade voting

Within each database the winning clade is chosen by a score-weighted vote: each domain contributes a length-normalized score (dom_score / model_len) to its clade, and the highest summed score wins:

$$\hat{c} = \arg\max_{c} \sum_{d \in D_c} \frac{\text{score}_d}{\text{len}_d}$$

Normalizing by model length controls for the tendency of longer profiles to accumulate higher raw scores, making domains comparable. Ties fall through to the existing mixture/completeness rules; Order and Superfamily are inherited from the winning clade.

This replaces TEsorter's raw domain-count plurality, which breaks ties by position in an internal collection rather than by any biological signal, causing sibling-clade swaps (Reina↔Tekay, Ale↔Alesia) to flip depending on search engine. Pass --compat-tesorter-voting to restore the original behaviour.

---

CLI reference

python3 src/pipeline.py <sequence> [options]

Flag	Default	Description
sequence	—	Input FASTA (TE library, or genome with --genome)
-d, --database	rexdb	Comma-separated database aliases or paths
--max-search	off	Search against all bundled databases
-o, --outdir	{input}.TEBinSorter	Output directory
--prefix	input basename	Output file prefix
-p, --processors	4	Processors
--facet	off	Facet pre-screen mode (AA databases only)
--bath	off	Frameshift-aware BATH engine (AA databases only)
--genome	off	Genome mode: domain-level annotation + GFF3
--win-size	1e6	Genome mode window size
--win-ovl	1e5	Genome mode window overlap
--emit-bath	off	Emit routed FASTA partitions for BATH to {outdir}/BATHwater/
--include-sine-so	off	Include the SINE_SO model in AnnoSINE searches
--compat-tesorter-voting	off	Raw domain-count clade vote (TEsorter behaviour)
--compat-tesorter-rounding	off	Round normalized scores to 2 dp before filtering (replicates a TEsorter bug)
--compat-tesorter-output	off	Emit combined .cls.tsv in TEsorter's 7-column format

---

Output files

File	Description
{prefix}.db	SQLite database with all hits, classifications and BLAST results
{prefix}.aa	Six-frame translated amino-acid sequences (indexed; HMMER path only)
{prefix}.{db}.cls.tsv	Per-database classifications (order, superfamily, clade, completeness)
{prefix}.cls.tsv	Combined classifications across databases + BLAST pass-2 (+ SecondaryHits)
{prefix}.{db}.classifications.tsv	Facet classifications with confidence tiers (--facet)
{prefix}.dom.gff3	Genome mode: classified TE protein-domain features
{prefix}.dom.faa / .dom.fna	Genome mode: domain sequences (AA for HMMER, nucleotide for BATH)
{prefix}.genome.summary.tsv	Genome mode: Order/Superfamily/Clade tallies
blast_pass2/	BLAST database and query chunks (temporary)
cut.fa	Genome mode: windowed genome (temporary)

---

Benchmarks

All runs on an ANVIL CPU node (Rosen Center for Advanced Computing, Purdue University; AMD EPYC 7763, 256 GB RAM), 16 threads, mean ± SD over 3 replicates. BATH in frameshift-aware mode (--fs).

Element mode

108,318 RepBase elements (LTR, TIR, LINE), each searched against its order-specific database (GyDB, REXdb-pnas, REXdb-line).

Pipeline	Engine	Time	Speedup vs TEsorter
TEsorter	HMMER (hmmscan)	2,932 s	1.0×
TEsorter2	pyHMMER	543 s	5.4×
TEsorter2	BATH	983 s	3.0×

Genome mode

Complete Oryza sativa genome (~375 Mb) against REXdb.

Pipeline	Engine	Time	Speedup vs TEsorter
TEsorter	HMMER (hmmscan)	2,127 s	1.0×
TEsorter2	pyHMMER	988 s	2.2×
TEsorter2	BATH	431 s	4.9×

BATH is 2.3× faster than pyHMMER here because it avoids the six-frame translation and coordinate back-mapping that dominate the HMMER path on long sequences.

---

Output files

File	Description
{prefix}.db	SQLite database with all hits, classifications and BLAST results
{prefix}.aa	Six-frame translated amino-acid sequences (indexed; HMMER path only)
{prefix}.{db}.cls.tsv	Per-database classifications (order, superfamily, clade, completeness)
{prefix}.cls.tsv	Combined classifications across databases + BLAST pass-2 (+ SecondaryHits)
{prefix}.{db}.classifications.tsv	Facet classifications with confidence tiers (--facet)
{prefix}.dom.gff3	Genome mode: classified TE protein-domain features
{prefix}.dom.faa / .dom.fna	Genome mode: domain sequences (AA for HMMER, nucleotide for BATH)
{prefix}.genome.summary.tsv	Genome mode: Order/Superfamily/Clade tallies
blast_pass2/	BLAST database and query chunks (temporary)
cut.fa	Genome mode: windowed genome (temporary)

---

TEsorter2 | BATH | 983 s | 3.0× |

---

TEsorter compatibility

src/tesorter_compat.py provides a drop-in CLI with TEsorter's original argument names and defaults (including count-based clade voting), for substituting TEsorter inside existing pipelines:

python3 src/tesorter_compat.py input.fasta -db rexdb -p 16 -pre out

Supported: -db/--hmm-database, --db-hmm, -st/--seq-type, -pre/--prefix, -p/--processors, -tmp/--tmp-dir, -cov/--min-coverage, -eval/--max-evalue, -prob/--min-probability, -score/--min-score, -dp2/--disable-pass2, -nolib/--no-library, -norc/--no-reverse, -nocln/--no-cleanup, plus --facet.

---

Architecture

Core

• pipeline.py — CLI and search orchestration
• search.py — HMM search engine with balanced parallelism
• sequence.py — FASTA ingestion (pyfastx) and six-frame translation (pyhmmer)
• hmm.py — HMM loading, alphabet detection, optimized profile construction
• bath_search.py — BATH engine: conversion, bathsearch --fs invocation, hit normalization
• genome.py — Genome mode: windowing, per-domain classification, overlap resolution, GFF3/summary
• results.py — SQLite persistence with pre-parsed columns (base_seq, strand, frame, domain_type)
• deconflict.py — numpy-based hit deconfliction and parameterized filtering

Facet classification

• decompose_hmm.py — standalone sub-HMM decomposition and splicing (tiered windows, configurable overlap; works for DNA and AA HMMs)
• model_graph.py — cross-model similarity graph, precomputed for bundled databases
• facet_classify.py — screen → verify → cross-family completion → legacy fallback
• cross_family.py — targeted search for missing domain families in classified frames

Classification and post-processing

• classifier.py — config-driven classification from domain hits; per-database domain remapping, overlap-aware deconfliction, order/superfamily/clade assignment
• blast_pass2.py — parallel chunked BLAST pass-2 with cross-database target pooling
• tesorter_output.py, emit.py, id_registry.py — output formatting and identifier bookkeeping

---

Extended methods

HMM facets

pyHMMER exposes an HMM's emission probabilities in Python. TEsorter2 uses them to locate conserved subregions that most influence HMMER's decision-making, and extracts each such region into a "facet": a complete, self-contained HMM whose emission and transition probabilities are cloned from its parent over the corresponding window.

Facets are sized at 96, 64, 48 or 32 amino acids, using the longest size the parent supports; models of ≤32 aa are used as-is. Windows are chosen to maximize the summed emission probability over the parent and may overlap by at most 33% of facet length.

For a parent HMM of length M and a query of length N, one hmmsearch is a dynamic-programming matrix of size M × N. A facet search replaces M with a facet size F (F < M), yielding an F × N slice with a proportional reduction in work. Three properties make this profitable:

• Short models are more decisive. A full-length model accumulates information about search effort across the whole sequence; a short model confirms or rejects a local region quickly.
• Facet scores predict full-length scores. A good facet hit almost always implies a good full-length hit, making facets a fast approximation of the final score.
• Facet sizes pack SIMD lanes. 96/64/32 aa are consumed in 16-aa bites by HMMER's internals, reducing low-level CPU waste relative to less divisible sizes. Facets do not reproduce their parent's domain detections exactly, so a full-length verification is still required for correctness. The acceleration comes from skipping verifications: the TE HMMs within each database are highly redundant (either wholly, as in single-type databases, or in subcollections, as in REXdb and GyDB), so most models in a cluster hit the same sequence at differing strengths and all but the best are discarded downstream anyway — yet detecting a weak hit costs exactly as much as detecting a strong one. TEsorter2 ships precomputed all-vs-all similarity graphs (obtained by searching each model's consensus against every other model in the database) that identify these clusters. Facet hits order the parent searches within each cluster from highest to lowest, and verification stops as soon as a high-scoring parent is confirmed — typically 1–2 parent searches per cluster per sequence instead of dozens.

In effect, the deconfliction that would otherwise happen after an exhaustive search is moved before it, at the cost of a cheap approximate screen. Post-filter agreement with the exhaustive search is 99.98% at hit level and 99.8% at family recall (4 misses of 2,090 on rice REXdb).

The facet generation code is intentionally a separate module and is reusable in other projects.

Parallel strategy

Default search. pyHMMER exposes two C-level parallelization schemes: queries, where each thread takes one HMM and searches it against all sequences, and targets, where one sequence is searched against models in parallel. queries is inherently more efficient unless there are few models.

HMM runtime scales roughly with M², so a single long model can dominate. AnnoSINE is the pathological case: SINE_SO (M≈4,100) accounts for ~71% of the model set's runtime, and under queries parallelism every other model finishes quickly while SINE_SO runs single-threaded for

10× longer than all the rest combined.

TEsorter2 precomputes each model's expected cost (M²) and bins them: small models (cost ≤ 75th percentile + 2×IQR for that database) run in queries mode; large models run in targets mode. Sequences are reused from the same in-memory object across both searches, so the split is essentially free, yielding near-full CPU utilization in the most efficient mode available for each model class.

Facet search. Staged: facets from all models are searched with --nobias exactly as in the default search; top facet hits per sequence are ordered by quality, noting each facet's parent; each sequence is searched against its best facet's parent at full length. If it verifies, no further searches run; otherwise the next-best facet hit for a different parent is tried. Most sequences verify on the first attempt; almost all within three. Because verification stops at a single best hit, a cross-family pass then searches each verified sequence against the top facet hit's parent for every other TE family, preserving both primary and secondary labelling sensitivity. Sequences still unclassified fall back to the legacy search.

The expensive part of HMMER is finding a match. Facets route each sequence only to the models likely to produce its best hits, skipping the redundant weak hits that would be discarded anyway. The leftovers sent to the legacy search are mostly genuine rejects with no TE signal, and are cheap: a REXdb legacy search spends ~92% of its runtime on the ~71k of ~880k protein frames that actually contain hits.

---