Add optional minimap2 pass-2 aligner; default behavior unchanged

README.md

@@ -1,3 +1,85 @@
+# TEsorter2_minimap2
+
+Fork of TEsorter2 that adds a [minimap2](https://github.com/lh3/minimap2)
+option for the pass-2 similarity search. By default, pass-2 behaves exactly
+like upstream TEsorter2 (`blastn`, 80-80-80 thresholds); opt into minimap2
+with `--pass2-aligner minimap2`. The minimap2 path runs a sensitivity-tuned
+flag set, then reduces the PAF to one row per query via
+`classify_ltr_paf_fast`, which enforces **identity, qcov, and tcov together**
+under the user-supplied I-C-L rule.
+
+## Additional runtime dependency
+
+`minimap2` binary must be on `$PATH` (only needed with
+`--pass2-aligner minimap2`). Install with conda:
+
+```
+mamba install -c bioconda minimap2
+```
+
+Everything else is unchanged from TEsorter2 (pyhmmer, pyfastx, numpy).
+
+## New / changed CLI options
+
+| Option | Default | Purpose |
+|---|---|---|
+| `--pass2-aligner {blast,minimap2}` | `blast` | Pass-2 engine. `blast` reproduces upstream TEsorter2's blastn pass-2; `minimap2` uses the PAF qcov+tcov path |
+| `--blast-task {megablast,dc-megablast}` | `megablast` | blastn `-task` for the blast engine. `dc-megablast` is slower but more sensitive to diverged matches |
+| `-dp2`, `--disable-pass2` | off | Skip pass-2 (HMM-only classification) |
+| `-rule`, `--pass2-rule I-C-L` | `80-80-80` | Pass-2 threshold. blast: pident, qcovs, and alignment-length filters. minimap2: I drives `--min-pid`, C drives both `--min-qcov` and `--min-tcov`, L is parsed for grammar compatibility but is not consumed by `classify_ltr_paf_fast` |
+| `--pass2-classified-fasta FASTA` | none | Optional FASTA of prior classifications to augment the pass-2 target pool. Headers must be shaped `>id#Order/Superfamily/Clade` |
+| `--minimap2-preset PRESET` | `asm20` | Passed through as `minimap2 -x` |
+| `--minimap2-extra STR` | empty | Additional flags appended to the minimap2 command line |
+
+## minimap2 invocation
+
+With `--pass2-aligner minimap2`, pass-2 runs (with target =
+previously-classified pool, query = HMM-unclassified):
+
+```
+minimap2 -x asm20 --rmq=no --no-long-join \
+    -k 10 -w 10 -r 500,20000 -g 500 -p 0.3 -N 100 -m 30 \
+    -t NCPU -K 1G --seed 11 --paf-no-hit \
+    -o pass2.paf  TARGET.fa  QUERY.fa
+```
+
+The PAF is then collapsed by `classify_ltr_paf_fast` to one row per query:
+
+```
+qname    pass/fail    pid    eff_qcov    eff_tcov    best_tname
+```
+
+A query is rescued (best target's classification inherited) iff the row
+reads `pass`, i.e. **pid ≥ I/100, eff_qcov ≥ C/100, eff_tcov ≥ C/100**.
+
+Benchmark at `70-70-70`: F1 ≈ 0.895, accuracy ≈ 0.943, precision ≈ 0.866,
+recall ≈ 0.926, MCC ≈ 0.857.
+
+## SQLite schema
+
+The `blast_hits` table stores one row per query (best target) with columns
+`qseqid, sseqid, pident, qcovs, tcovs, passes_rule, classified_by`. Indexes
+on `qseqid` and `sseqid` are still built by `results.finalize_db`.
+
+## What changed vs stock TEsorter2
+
+- `tesorter2/blast_pass2.py` — internals swapped from `blastn`+`multiprocessing.Pool`
+  to one `minimap2` call followed by `classify_ltr_paf_fast.process_paf`.
+- `tesorter2/minimap.py` — minimap2 wrapper with the sensitivity-tuned pass-2 flag
+  set (no `-c`; relies on PAF + `dv:f` only).
+- `tesorter2/classify_ltr_paf_fast.py` — PAF → TSV reducer (one row per query, with
+  pass/fail under the I-C-L rule).
+- `tesorter2/pass2_external.py` — shared with mmseqs port. Helpers for
+  `--pass2-classified-fasta`.
+- `tesorter2/pipeline.py` / `tesorter2/tesorter_compat.py` — wire the five new CLI args
+  through the pass-2 call.
+
+---
+
+Upstream TEsorter2 README follows.
+
+---
+
 # TEsorter2
 
 Fast, divergence-robust classification of transposable elements.
@@ -14,6 +96,8 @@ keeps its classification semantics while introducing three major improvements:
    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.
 
+Note: I should switch from default task (megablast) to dc-megablast.
+
 ## Installation
 
 ### conda (recommended)

tesorter2/blast_backend.py

@@ -0,0 +1,259 @@
+"""
+blast_backend.py — TEsorter2 master's blastn pass-2 logic, lifted verbatim so
+the `--pass2-aligner blast` path runs and post-processes identically to
+TEsorter2 master.
+
+The orchestrator `run_pass2_blast` is called by blast_pass2.blast_pass2() AFTER
+the shared classified/unclassified split and the optional pass2_external merge,
+so the external-pool augmentation (`--pass2-classified-fasta`) applies to both
+aligner backends. Only the alignment + parse + classify is master-specific here.
+
+master functions (make_blast_db, run_blast_chunk, parse_blast_output,
+store_blast_hits, classify_from_blast) are copied from
+origin/master:src/blast_pass2.py; run_blast_chunk additionally takes a
+blast_task parameter (--blast-task). chunk_fasta + run_pass2_blast are new.
+"""
+
+import logging
+import multiprocessing
+import os
+import subprocess
+import time
+
+import pyfastx
+
+log = logging.getLogger(__name__)
+
+
+# ---- lifted verbatim from origin/master:src/blast_pass2.py ----
+
+def make_blast_db(db_fasta, seq_type="nucl"):
+    """Run makeblastdb."""
+    dbtype = seq_type
+    cmd = f"makeblastdb -in {db_fasta} -dbtype {dbtype} -out {db_fasta}"
+    result = subprocess.run(cmd, shell=True, capture_output=True, text=True)
+    if result.returncode != 0:
+        log.error(f"makeblastdb failed: {result.stderr}")
+        raise RuntimeError(f"makeblastdb failed: {result.stderr}")
+    log.info(f"  BLAST database built: {db_fasta}")
+
+
+def run_blast_chunk(query_chunk, db_fasta, output, seq_type="nucl", ncpu=1,
+                    blast_task="megablast"):
+    """Run BLAST on one query chunk."""
+    app = "blastn" if seq_type == "nucl" else "blastp"
+    outfmt = ("6 qseqid sseqid pident length mismatch gapopen qstart qend "
+              "sstart send evalue bitscore qlen slen qcovs qcovhsp sstrand")
+    # The -task value is blastn-only; blastp would reject it, so gate it on
+    # the blastn branch.
+    task = f" -task {blast_task}" if app == "blastn" else ""
+    cmd = (f"{app}{task} -query {query_chunk} -db {db_fasta} -out {output} "
+           f"-outfmt '{outfmt}' -num_threads {ncpu}")
+    result = subprocess.run(cmd, shell=True, capture_output=True, text=True)
+    if result.returncode != 0:
+        log.warning(f"BLAST chunk failed: {result.stderr[:200]}")
+    return output
+
+
+def parse_blast_output(blast_out):
+    """Parse BLAST outfmt 6 into hit dicts."""
+    fields = ["qseqid", "sseqid", "pident", "length", "mismatch", "gapopen",
+              "qstart", "qend", "sstart", "send", "evalue", "bitscore",
+              "qlen", "slen", "qcovs", "qcovhsp", "sstrand"]
+    types = [str, str, float, int, int, int, int, int, int, int,
+             float, float, int, int, float, float, str]
+
+    hits = []
+    if not os.path.exists(blast_out):
+        return hits
+
+    with open(blast_out) as f:
+        for line in f:
+            vals = line.strip().split("\t")
+            if len(vals) < len(fields):
+                continue
+            hit = {}
+            for field, typ, val in zip(fields, types, vals):
+                hit[field] = typ(val)
+            hits.append(hit)
+
+    return hits
+
+
+def store_blast_hits(conn, hits, db_seq_to_dbs):
+    """Store BLAST hits in SQLite."""
+    conn.execute("""
+        CREATE TABLE IF NOT EXISTS blast_hits (
+            qseqid      TEXT NOT NULL,
+            sseqid      TEXT NOT NULL,
+            pident      REAL NOT NULL,
+            length      INTEGER NOT NULL,
+            evalue      REAL NOT NULL,
+            bitscore    REAL NOT NULL,
+            qlen        INTEGER NOT NULL,
+            slen        INTEGER NOT NULL,
+            qcovs       REAL NOT NULL,
+            classified_by TEXT NOT NULL
+        )
+    """)
+
+    rows = []
+    for h in hits:
+        dbs = db_seq_to_dbs.get(h["sseqid"], set())
+        classified_by = ",".join(sorted(dbs)) if dbs else "unknown"
+        rows.append((
+            h["qseqid"], h["sseqid"], h["pident"], h["length"],
+            h["evalue"], h["bitscore"], h["qlen"], h["slen"],
+            h["qcovs"], classified_by,
+        ))
+
+    conn.executemany(
+        "INSERT INTO blast_hits VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?)",
+        rows,
+    )
+    conn.commit()
+    log.info(f"  Stored {len(rows)} BLAST hits")
+
+
+def classify_from_blast(conn, classifications, database=None,
+                        min_identity=80, min_coverage=80, min_length=80):
+    """Classify unclassified sequences from BLAST hits (master logic)."""
+    tables = {r[0] for r in conn.execute(
+        "SELECT name FROM sqlite_master WHERE type='table'").fetchall()}
+    if "blast_hits" not in tables:
+        return []
+
+    where = "WHERE pident >= ? AND qcovs >= ? AND length >= ?"
+    params = [min_identity, min_coverage, min_length]
+
+    if database:
+        where += " AND classified_by LIKE ?"
+        params.append(f"%{database}%")
+
+    rows = conn.execute(f"""
+        SELECT qseqid, sseqid, pident, qcovs, length, bitscore
+        FROM blast_hits
+        {where}
+        ORDER BY bitscore DESC
+    """, params).fetchall()
+
+    classified_set = set(classifications.keys())
+
+    best = {}
+    for qid, sid, pident, qcovs, length, bitscore in rows:
+        if qid in classified_set:
+            continue
+        if qid not in best:
+            best[qid] = (sid, pident, qcovs, length, bitscore)
+
+    new_classifications = []
+    no_source = 0
+    for qid, (sid, pident, qcovs, length, bitscore) in best.items():
+        if sid in classifications:
+            source = classifications[sid]
+            new_classifications.append({
+                "id": qid,
+                "order": source["order"],
+                "superfamily": source["superfamily"],
+                "clade": "unknown",
+                "complete": "none",
+                "strand": "?",
+                "domains": "none",
+                "blast_source": sid,
+                "blast_pident": pident,
+                "blast_qcovs": qcovs,
+                "blast_bitscore": bitscore,
+            })
+        else:
+            no_source += 1
+
+    if no_source:
+        log.info(f"    {no_source} BLAST hits to unclassified targets (skipped)")
+
+    log.info(f"  BLAST pass-2: {len(new_classifications)} sequences classified "
+             f"(from {len(best)} hits passing filters)")
+    return new_classifications
+
+
+# ---- new: chunking + orchestration ----
+
+def chunk_fasta(qry_fasta, n_chunks, outdir):
+    """Bin-pack sequences from qry_fasta into n_chunks files by total length.
+
+    Mirrors master's split bin-packing, but operates on an already-written
+    unclassified-query FASTA (the shared split in blast_pass2 produced it).
+    Returns the list of non-empty chunk paths. Returns [] for a missing or
+    empty input (pyfastx raises on empty files).
+    """
+    os.makedirs(outdir, exist_ok=True)
+    if not os.path.exists(qry_fasta) or os.path.getsize(qry_fasta) == 0:
+        return []
+    chunk_paths = [os.path.join(outdir, f"blast_query_{i}.fa")
+                   for i in range(max(1, n_chunks))]
+    handles = []
+    try:
+        handles = [open(p, "w") for p in chunk_paths]
+        lengths = [0] * len(chunk_paths)
+        fa = pyfastx.Fasta(qry_fasta, build_index=True)
+        for rec in fa:
+            i = lengths.index(min(lengths))
+            handles[i].write(f">{rec.name}\n{rec.seq}\n")
+            lengths[i] += len(rec.seq)
+    finally:
+        for h in handles:
+            h.close()
+    return [p for p in chunk_paths if os.path.getsize(p) > 0]
+
+
+def run_pass2_blast(qry_fasta, db_fasta, conn, classifications, db_seq_to_dbs,
+                    n_processors, min_identity, min_coverage, min_length, work,
+                    blast_task="megablast"):
+    """blastn pass-2 over an already-prepared (db_fasta, qry_fasta) pair.
+
+    Reproduces TEsorter2 master's: makeblastdb -> chunked parallel blastn ->
+    outfmt6 parse -> SQLite -> classify_from_blast (qcovs+length filter, best by
+    bitscore, clade=unknown). The I/C/L thresholds come from the run's -rule.
+    """
+    t0 = time.time()
+    os.makedirs(work, exist_ok=True)
+
+    if not os.path.exists(db_fasta) or os.path.getsize(db_fasta) == 0:
+        log.info("  pass-2 target FASTA is empty or missing; skipping blastn")
+        return []
+
+    make_blast_db(db_fasta, seq_type="nucl")
+
+    query_chunks = chunk_fasta(qry_fasta, n_processors, work)
+    if not query_chunks:
+        log.info("  No unclassified sequences to search")
+        return []
+
+    log.info(f"  Running {len(query_chunks)} BLAST processes")
+    t1 = time.time()
+    blast_outputs = []
+    args_list = []
+    for chunk in query_chunks:
+        out = chunk + ".blastout"
+        blast_outputs.append(out)
+        args_list.append((chunk, db_fasta, out, "nucl", 1, blast_task))
+    with multiprocessing.Pool(len(query_chunks)) as pool:
+        pool.starmap(run_blast_chunk, args_list)
+    t2 = time.time()
+    log.info(f"  BLAST search: {t2 - t1:.1f}s")
+
+    all_hits = []
+    for blast_out in blast_outputs:
+        all_hits.extend(parse_blast_output(blast_out))
+    log.info(f"  {len(all_hits)} total BLAST hits")
+
+    if all_hits:
+        store_blast_hits(conn, all_hits, db_seq_to_dbs)
+
+    new_cls = classify_from_blast(
+        conn, classifications,
+        min_identity=min_identity,
+        min_coverage=min_coverage,
+        min_length=min_length,
+    )
+    log.info(f"  BLAST pass-2 total: {time.time() - t0:.1f}s")
+    return new_cls