diff --git a/README.md b/README.md index 022fdf4..ab5ef26 100644 --- a/README.md +++ b/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) diff --git a/tesorter2/blast_backend.py b/tesorter2/blast_backend.py new file mode 100644 index 0000000..0b465d5 --- /dev/null +++ b/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 diff --git a/tesorter2/blast_pass2.py b/tesorter2/blast_pass2.py index 662cf12..77108f9 100644 --- a/tesorter2/blast_pass2.py +++ b/tesorter2/blast_pass2.py @@ -1,369 +1,299 @@ """ -blast_pass2.py — BLAST-based pass-2 classification for HMM-unclassified sequences. - -Sequences not classified by HMM search are BLASTed against classified sequences. -If a strong similarity match exists (80-80-80 rule by default), the unclassified -sequence inherits the classification of its best BLAST hit. - -Key design: - - Cross-database: one BLAST search against all classified sequences from all databases - - Per-database reconstruction via filtering on classified_by - - Parallel chunked BLAST with greedy bin-packing by sequence length - - All results stored in SQLite for post-hoc threshold adjustment +blast_pass2.py — minimap2-based pass-2 classification for HMM-unclassified +sequences. + +Pipeline: + 1. minimap2 (sensitivity-tuned flags from `minimap.run_minimap2`) writes a + PAF for the unclassified queries against the classified-pool target. + 2. `classify_ltr_paf_fast.process_paf` reduces the PAF to one row per query: + qname pass/fail pid eff_qcov eff_tcov best_tname + under the rule `--min-pid I --min-qcov C --min-tcov C` derived from + `--pass2-rule I-C-L`. (Per benchmarking, 70-70-70 is recommended for + this minimap2 path; the CLI default is 80-80-80 to mirror upstream's + blastn pass-2. The L value is parsed for backwards-compat with the + I-C-L grammar but is not consumed downstream.) + 3. Each `pass` row inherits the target's order/superfamily/clade. + +The SQLite `blast_hits` table is preserved for post-run introspection but the +schema has been simplified to the columns classify_ltr_paf_fast emits. """ import logging import os -import subprocess import tempfile import time from collections import defaultdict import pyfastx +from . import minimap +from . import pass2_external +from .classify_ltr_paf_fast import process_paf, format_row + log = logging.getLogger(__name__) def _get_classified_ids(conn): - """Get classified sequence IDs per database from classifier results. - - Returns: - classified: {base_seq: set(databases)} — which databases classified each seq - """ + """Get classified sequence IDs per database from classifier results.""" classified = defaultdict(set) - - # Check which tables exist tables = {r[0] for r in conn.execute( "SELECT name FROM sqlite_master WHERE type='table'").fetchall()} - if "legacy_hits" in tables: for row in conn.execute( "SELECT DISTINCT base_seq, database FROM legacy_hits" ): classified[row[0]].add(row[1]) - return dict(classified) def split_classified_unclassified(input_fasta, classified_ids, outdir, - n_chunks=4, seq_type="nucl"): - """Split input into classified (BLAST db) and chunked unclassified (queries). - - One pass through the input. Unclassified sequences are bin-packed into - n_chunks files by total sequence length for even BLAST parallelism. - - Args: - input_fasta: path to input FASTA - classified_ids: {seq_name: set(databases)} - outdir: directory for output files - n_chunks: number of query chunks - seq_type: "nucl" or "prot" (determines which sequences to write) - - Returns: - db_fasta: path to classified sequences FASTA (BLAST database) - query_chunks: list of paths to unclassified sequence chunks - db_seq_to_dbs: {seq_name: set(databases)} for DB sequences + seq_type="nucl"): + """Split input into classified-pool FASTA (target) and unclassified-query + FASTA. Nucleotide pools are uppercased and stripped of non-ATCG to keep + inputs clean (minimap2 itself tolerates ambiguous bases). """ os.makedirs(outdir, exist_ok=True) - db_fasta = os.path.join(outdir, "blast_db.fa") - chunk_paths = [os.path.join(outdir, f"blast_query_{i}.fa") for i in range(n_chunks)] - - # Open all handles - db_handle = open(db_fasta, "w") - chunk_handles = [open(p, "w") for p in chunk_paths] - chunk_lengths = [0] * n_chunks + qry_fasta = os.path.join(outdir, "blast_query.fa") + nucl = (seq_type == "nucl") fa = pyfastx.Fasta(input_fasta, build_index=True) + n_classified = 0 n_unclassified = 0 db_seq_to_dbs = {} - for rec in fa: - name = rec.name - seq = str(rec.seq) + with open(db_fasta, "w") as dbh, open(qry_fasta, "w") as qh: + for rec in fa: + name = rec.name + seq = str(rec.seq) + if nucl: + seq = "".join(c for c in seq.upper() if c in "ATCG") + if name in classified_ids: + dbh.write(f">{name}\n{seq}\n") + db_seq_to_dbs[name] = classified_ids[name] + n_classified += 1 + else: + qh.write(f">{name}\n{seq}\n") + n_unclassified += 1 + + log.info(f" Split: {n_classified} classified (DB), " + f"{n_unclassified} unclassified (query)") + return db_fasta, qry_fasta, db_seq_to_dbs - if name in classified_ids: - db_handle.write(f">{name}\n{seq}\n") - db_seq_to_dbs[name] = classified_ids[name] - n_classified += 1 - else: - # Bin-pack into lightest chunk - min_idx = chunk_lengths.index(min(chunk_lengths)) - chunk_handles[min_idx].write(f">{name}\n{seq}\n") - chunk_lengths[min_idx] += len(seq) - n_unclassified += 1 - db_handle.close() - for h in chunk_handles: - h.close() +def run_alignment(query_fa, target_fa, paf_out, ncpu=4, + preset="asm20", extra=""): + """Run minimap2 with the sensitivity-tuned pass-2 flag set.""" + minimap.run_minimap2( + query_fa=query_fa, target_fa=target_fa, paf_out=paf_out, + ncpu=ncpu, preset=preset, extra=extra, + ) + return paf_out - # Remove empty chunks - query_chunks = [p for p in chunk_paths if os.path.getsize(p) > 0] - log.info(f" Split: {n_classified} classified (DB), " - f"{n_unclassified} unclassified -> {len(query_chunks)} chunks") - - return db_fasta, query_chunks, db_seq_to_dbs - - -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): - """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" - cmd = (f"{app} -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. - - Adds classified_by field indicating which databases classified - each target sequence. +def classify_paf_to_tsv(paf_path, tsv_out, min_pid, min_qcov, min_tcov): + """Run classify_ltr_paf_fast over the PAF, write the TSV alongside, + and return the parsed rows: [(qname, pass_str, pid, qcov, tcov, best_tname), ...]. """ - # Indexes on blast_hits are built by results.finalize_db at end of run. + if not os.path.exists(paf_path) or os.path.getsize(paf_path) == 0: + return [] + + with open(paf_path) as fh: + rows = process_paf( + fh, min_pid=min_pid, min_qcov=min_qcov, min_tcov=min_tcov, + fill_colinear=False, + verbose=False, + ) + + with open(tsv_out, "w") as fh: + for row in rows: + fh.write(format_row(row) + "\n") + + n_pass = sum(1 for r in rows if r[1] == "pass") + log.info(f" classify_ltr_paf_fast: {n_pass}/{len(rows)} queries pass " + f"(min_pid={min_pid:.3f} min_qcov={min_qcov:.3f} " + f"min_tcov={min_tcov:.3f}) -> {tsv_out}") + return rows + + +def store_blast_hits(conn, tsv_rows, db_seq_to_dbs): + """Store classify_ltr_paf_fast rows in SQLite. One row per query.""" 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, + qseqid TEXT NOT NULL, + sseqid TEXT NOT NULL, + pident REAL NOT NULL, + qcovs REAL NOT NULL, + tcovs REAL NOT NULL, + passes_rule INTEGER NOT NULL, classified_by TEXT NOT NULL ) """) - rows = [] - for h in hits: - dbs = db_seq_to_dbs.get(h["sseqid"], set()) + for qname, pass_str, pid, qcov, tcov, tname in tsv_rows: + dbs = db_seq_to_dbs.get(tname, 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, + qname, tname, + pid * 100.0, qcov * 100.0, tcov * 100.0, + 1 if pass_str == "pass" else 0, + classified_by, )) - conn.executemany( - "INSERT INTO blast_hits VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?)", - rows, + "INSERT INTO blast_hits VALUES (?, ?, ?, ?, ?, ?, ?)", rows, ) conn.commit() - log.info(f" Stored {len(rows)} BLAST hits") + log.info(f" Stored {len(rows)} minimap2 best-hit rows") -def classify_from_blast(conn, classifications, database=None, - min_identity=80, min_coverage=80, min_length=80): - """Classify unclassified sequences from BLAST hits. - - Args: - conn: sqlite3 connection with blast_hits table - classifications: dict of {seq_id: {order, superfamily, ...}} from - classifier.classify_sequences() - database: if set, only accept BLAST targets classified by this database - min_identity: minimum percent identity - min_coverage: minimum query coverage - min_length: minimum alignment length - - Returns: - list of classification dicts for newly classified sequences +def classify_from_blast(tsv_rows, classifications): + """Inherit order/superfamily/clade for queries whose best target passed + the rule. `classifications` maps target_id -> classification dict. """ - # Check table exists - tables = {r[0] for r in conn.execute( - "SELECT name FROM sqlite_master WHERE type='table'").fetchall()} - if "blast_hits" not in tables: - return [] - - # Build filter - 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}%") - - # Best hit per query by bitscore - rows = conn.execute(f""" - SELECT qseqid, sseqid, pident, qcovs, length, bitscore - FROM blast_hits - {where} - ORDER BY bitscore DESC - """, params).fetchall() - - # Classified ID set for quick lookup classified_set = set(classifications.keys()) - - best = {} - for qid, sid, pident, qcovs, length, bitscore in rows: - if qid in classified_set: - continue # already classified by HMM, skip - if qid not in best: - best[qid] = (sid, pident, qcovs, length, bitscore) - - # Inherit classification from best hit's target 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: + n_pass = 0 + for qname, pass_str, pid, qcov, tcov, tname in tsv_rows: + if pass_str != "pass": + continue + n_pass += 1 + if qname in classified_set: + continue + if tname not in classifications: no_source += 1 + continue + source = classifications[tname] + new_classifications.append({ + "id": qname, + "order": source["order"], + "superfamily": source["superfamily"], + "clade": source.get("clade", "unknown"), + "complete": "none", + "strand": "?", + "domains": "none", + "blast_source": tname, + "blast_pident": pid * 100.0, + "blast_qcovs": qcov * 100.0, + "blast_tcovs": tcov * 100.0, + "blast_bitscore": 0.0, + }) 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)") + log.info(f" {no_source} pass-2 hits to unclassified targets (skipped)") + log.info(f" pass-2: {len(new_classifications)} sequences classified " + f"(from {n_pass} rule-passing queries)") return new_classifications def blast_pass2(input_fasta, conn, hmm_classifications=None, seq_type="nucl", n_processors=4, min_identity=80, min_coverage=80, min_length=80, - outdir=None): - """Full BLAST pass-2 pipeline. + outdir=None, + pass2_classified_fasta=None, + preset="asm20", minimap2_extra="", + aligner="blast", blast_task="megablast"): + """Pass-2 similarity search (blastn by default, minimap2 opt-in). Args: - input_fasta: path to input FASTA - conn: sqlite3 connection with HMM results - hmm_classifications: dict of {seq_id: {order, superfamily, clade, ...}} - from classifier.classify_sequences(). Targets inherit - classification from their best BLAST match. - seq_type: "nucl" or "prot" - n_processors: number of parallel BLAST processes - min_identity: filter threshold - min_coverage: filter threshold - min_length: filter threshold - outdir: output directory (default: tempdir) - - Returns: - list of new classification dicts + min_identity: I from --pass2-rule I-C-L (percent, e.g. 80) + min_coverage: C from --pass2-rule I-C-L (percent; blast applies it to + qcovs, minimap2 to qcov AND tcov) + min_length: L from --pass2-rule I-C-L (blast: minimum alignment + length; minimap2: parsed for backwards-compat with the + I-C-L grammar, not consumed by classify_ltr_paf_fast) """ t0 = time.time() + if aligner == "minimap2": + minimap.check_minimap2() + + if seq_type != "nucl": + log.warning("minimap2 pass-2 only supports nucleotide sequences; " + f"seq_type={seq_type!r} will be treated as nucl") + seq_type = "nucl" - # Get classified IDs from database classified_ids = _get_classified_ids(conn) if not classified_ids: - log.info(" No classified sequences for BLAST pass-2") + log.info(" No classified sequences for pass-2") return [] - log.info(f" BLAST pass-2: {len(classified_ids)} classified sequences as targets") + log.info(f" pass-2 ({aligner}): {len(classified_ids)} classified sequences as targets") - # Split if outdir is None: - outdir = tempfile.mkdtemp(prefix="tesorter2_blast_") + outdir = tempfile.mkdtemp(prefix="tesorter2_minimap2_") + work = os.path.join(outdir, "minimap2_pass2") - blast_dir = os.path.join(outdir, "blast_pass2") - db_fasta, query_chunks, db_seq_to_dbs = split_classified_unclassified( - input_fasta, classified_ids, blast_dir, n_chunks=n_processors, - seq_type=seq_type) + db_fasta, qry_fasta, db_seq_to_dbs = split_classified_unclassified( + input_fasta, classified_ids, work, seq_type=seq_type) - if not query_chunks: + if os.path.getsize(qry_fasta) == 0: log.info(" No unclassified sequences to search") return [] - # Build BLAST database - make_blast_db(db_fasta, seq_type=seq_type) + if hmm_classifications is None: + hmm_classifications = {} - # Run parallel BLAST - log.info(f" Running {len(query_chunks)} BLAST processes") - t1 = time.time() + if pass2_classified_fasta: + updated = pass2_external.update_classified_fasta_headers( + pass2_classified_fasta, hmm_classifications, work + ) + src = updated or pass2_classified_fasta + pass2_external.extend_hmm_classifications_from_fasta( + hmm_classifications, src, db_seq_to_dbs + ) + merged_db = os.path.join(work, "pass2_db_merged.fa") + pass2_external.merge_classified_fastas( + merged_db, db_fasta, src, clean_nucl=True + ) + db_fasta = merged_db - import multiprocessing - blast_outputs = [] - args_list = [] - for i, chunk in enumerate(query_chunks): - out = chunk + ".blastout" - blast_outputs.append(out) - args_list.append((chunk, db_fasta, out, seq_type, 1)) + if os.path.getsize(db_fasta) == 0: + log.info(" pass-2 target FASTA is empty; skipping pass-2") + return [] - with multiprocessing.Pool(n_processors) as pool: - pool.starmap(run_blast_chunk, args_list) + if aligner == "blast": + from . import blast_backend + new_cls = blast_backend.run_pass2_blast( + qry_fasta=qry_fasta, db_fasta=db_fasta, conn=conn, + classifications=hmm_classifications, + db_seq_to_dbs=db_seq_to_dbs, + n_processors=n_processors, + min_identity=min_identity, + min_coverage=min_coverage, + min_length=min_length, + work=work, + blast_task=blast_task, + ) + log.info(f" blast pass-2 total: {time.time() - t0:.1f}s") + return new_cls + + paf_out = os.path.join(work, "pass2.paf") + tsv_out = os.path.join(work, "pass2.tsv") + log.info(f" Running minimap2 -x {preset} with {n_processors} threads") + t1 = time.time() + run_alignment( + query_fa=qry_fasta, target_fa=db_fasta, paf_out=paf_out, + ncpu=n_processors, preset=preset, extra=minimap2_extra, + ) t2 = time.time() - log.info(f" BLAST search: {t2 - t1:.1f}s") - - # Parse and store - 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) - - # Build HMM classifications lookup for inheritance - hmm_cls = {} - if hmm_classifications is not None: - hmm_cls = hmm_classifications - log.info(f" {len(hmm_cls)} HMM classifications available for inheritance") - - # Classify - new_cls = classify_from_blast( - conn, hmm_cls, - min_identity=min_identity, - min_coverage=min_coverage, - min_length=min_length, + log.info(f" minimap2 alignment: {t2 - t1:.1f}s") + + min_pid = min_identity / 100.0 + min_cov = min_coverage / 100.0 + tsv_rows = classify_paf_to_tsv( + paf_out, tsv_out, + min_pid=min_pid, min_qcov=min_cov, min_tcov=min_cov, ) - t3 = time.time() - log.info(f" BLAST pass-2 total: {t3 - t0:.1f}s") + if tsv_rows: + store_blast_hits(conn, tsv_rows, db_seq_to_dbs) + log.info(f" {len(hmm_classifications)} HMM classifications available for inheritance") + + new_cls = classify_from_blast(tsv_rows, hmm_classifications) + + t3 = time.time() + log.info(f" minimap2 pass-2 total: {t3 - t0:.1f}s") return new_cls diff --git a/tesorter2/classify_ltr_paf_fast.py b/tesorter2/classify_ltr_paf_fast.py new file mode 100755 index 0000000..ae9eeee --- /dev/null +++ b/tesorter2/classify_ltr_paf_fast.py @@ -0,0 +1,413 @@ +#!/usr/bin/env python3 +"""Classify putative LTR-RTs from a minimap2 PAF -- CIGAR-free version. + +This is a faster, memory-leaner sibling of classify_ltr_paf.py that uses ONLY +the standard PAF columns plus the dv:f tag (approximate per-base sequence +divergence). It does NOT walk the cg:Z CIGAR. The intent is to let you run +minimap2 WITHOUT -c, saving substantial runtime on large inputs: + + minimap2 -k15 -w5 -A2 -B3 -r1k,10k -s30 -m30 -N50 -p0.1 \\ + target.fa query.fa > out.paf # NOTE: no -c, no cg:Z, dv:f tag + +Approximation tradeoff: q- and t-intervals for overlap dedup are taken as +[qstart, qend] and [tstart, tend], i.e. the alignment SPAN, which lumps +gap-induced bases (insertions and deletions) in with the matches+mismatches. +For LTR-RTs at >=70% identity, gap content is typically 2-5%, so eff_qcov / +eff_tcov are inflated by that small amount. dv:f is used in place of de:f +(de:f is gap-compressed; dv:f is approximate per-base divergence -- close +enough at our 70% threshold). + +If both dv:f and de:f are present (e.g. the user kept -c for some reason), +de:f is used. Otherwise dv:f. If neither is present, the line is skipped. + +OUTPUT (TSV; --header to add a header): + qname pass pid eff_qcov eff_tcov best_tname + +PARSIMONY RULE for pid (same as the CIGAR-aware sibling): + 1. Drop alignments fully encompassed by another on q (encompassing aln has + more statistical power). + 2. Sort survivors by dv:f DESCENDING (densest first). + 3. Greedy allocation, densest-first: each alignment's mutations are first + absorbed into already-claimed (denser) overlap; leftovers spill into + its unique region. + 4. pid = 1 - sum(unique_mutations) / sum(unique_q_spans). + +Per-alignment "mutations" weight: dv * qspan (where qspan = qend - qstart). + +Best-target per query: passes-rule first, then joint_score = pid*qcov*tcov, +then tname alphabetical. +""" + +import argparse +import sys +from collections import defaultdict + + +def merge_intervals(ivs): + """Merge overlapping/adjacent intervals; return sorted, disjoint list.""" + if not ivs: + return [] + ivs = sorted(ivs) + out = [list(ivs[0])] + for s, e in ivs[1:]: + if s <= out[-1][1]: + if e > out[-1][1]: + out[-1][1] = e + else: + out.append([s, e]) + return [(s, e) for s, e in out] + + +def interval_total(ivs): + return sum(e - s for s, e in ivs) + + +def interval_difference(A, B): + """Return A \\ B: parts of A not covered by B.""" + A = merge_intervals(A) + B = merge_intervals(B) + out = [] + bi = 0 + for s, e in A: + cur = s + while bi < len(B) and B[bi][1] <= cur: + bi += 1 + j = bi + while cur < e and j < len(B): + bs, be = B[j] + if bs >= e: + break + if bs > cur: + out.append((cur, min(bs, e))) + cur = max(cur, be) + j += 1 + if cur < e: + out.append((cur, e)) + return out + + +def interval_intersection_length(A, B): + A = merge_intervals(A) + B = merge_intervals(B) + i = j = 0 + total = 0 + while i < len(A) and j < len(B): + s = max(A[i][0], B[j][0]) + e = min(A[i][1], B[j][1]) + if s < e: + total += e - s + if A[i][1] < B[j][1]: + i += 1 + else: + j += 1 + return total + + +def is_fully_encompassed(small_iv, big_iv): + """True iff every interval in small_iv is fully inside some interval in big_iv.""" + big_iv = merge_intervals(big_iv) + for s, e in small_iv: + contained = False + for bs, be in big_iv: + if bs <= s and e <= be: + contained = True + break + if not contained: + return False + return True + + +def parse_paf_line(line, lineno, strict=False): + """Return alignment record dict, or None if malformed and strict=False. + + Required: 12 standard columns + a dv:f or de:f tag. + """ + fields = line.rstrip("\n").split("\t") + if len(fields) < 12: + if strict: + raise ValueError("line %d: only %d fields (need >=12)" % (lineno, len(fields))) + return None + dv = None + de = None + for tag in fields[12:]: + if tag.startswith("dv:f:"): + try: + dv = float(tag[5:]) + except ValueError: + pass + elif tag.startswith("de:f:"): + try: + de = float(tag[5:]) + except ValueError: + pass + # Prefer the gap-compressed de:f if present (more accurate); fall back to dv:f. + div = de if de is not None else dv + if div is None: + if strict: + raise ValueError("line %d: neither dv:f nor de:f tag present" % lineno) + return None + try: + return { + "qname": fields[0], + "qlen": int(fields[1]), + "qstart": int(fields[2]), + "qend": int(fields[3]), + "strand": fields[4], + "tname": fields[5], + "tlen": int(fields[6]), + "tstart": int(fields[7]), + "tend": int(fields[8]), + "div": div, + } + except (ValueError, IndexError): + if strict: + raise + return None + + +def chained_intervals(alignments, max_gap=5000, gap_tol=0.20): + """Group alignments into colinear chains; return per-chain (q,t) envelopes. + + Two adjacent alignments (sorted by qstart) join the same chain iff: + 1. Same strand. + 2. t-order consistent with strand (ascending t for +, descending t for -). + 3. Inner gap on each axis <= max_gap (bp). + 4. |q_gap - t_gap| / max(q_gap, t_gap) <= gap_tol (synchronized indel). + Otherwise a new chain begins. + + Returns ([(qmin, qmax), ...], [(tmin, tmax), ...]) -- one (q,t) envelope per + chain. Multiple chains may overlap on q or t; the caller is responsible for + merging across chains. + """ + if not alignments: + return [], [] + alns = sorted(alignments, key=lambda a: (a["qstart"], a["qend"])) + chains = [[alns[0]]] + for nxt in alns[1:]: + cur = chains[-1][-1] + join = cur["strand"] == nxt["strand"] + if join: + if cur["strand"] == "+": + t_order_ok = nxt["tstart"] >= cur["tstart"] + t_gap = max(0, nxt["tstart"] - cur["tend"]) + else: + t_order_ok = nxt["tstart"] <= cur["tstart"] + t_gap = max(0, cur["tstart"] - nxt["tend"]) + join = t_order_ok + if join: + q_gap = max(0, nxt["qstart"] - cur["qend"]) + big = max(q_gap, t_gap) + if big > max_gap: + join = False + elif big > 0 and abs(q_gap - t_gap) / big > gap_tol: + join = False + if join: + chains[-1].append(nxt) + else: + chains.append([nxt]) + q_ivs, t_ivs = [], [] + for chain in chains: + q_ivs.append((min(a["qstart"] for a in chain), + max(a["qend"] for a in chain))) + t_ivs.append((min(a["tstart"] for a in chain), + max(a["tend"] for a in chain))) + return q_ivs, t_ivs + + +def compute_pair_metrics(alignments, qlen, tlen, + fill_colinear=False, max_gap=5000, gap_tol=0.20): + """Compute (pid, eff_qcov, eff_tcov) for one (q,t) pair, span-based. + + `alignments` is a non-empty list of dicts: {qstart, qend, tstart, tend, div}. + If fill_colinear is True, colinear-chain inner gaps on q and t are bridged + in the qcov/tcov computation (pid logic is unchanged either way). + """ + # Build per-alignment span intervals (no CIGAR -- one interval each). + for a in alignments: + a["q_iv"] = [(a["qstart"], a["qend"])] + a["t_iv"] = [(a["tstart"], a["tend"])] + a["qspan"] = a["qend"] - a["qstart"] + + # 1. Drop alignments whose q-span is fully encompassed by a STRICTLY longer one. + survivors = [] + for i, A in enumerate(alignments): + encompassed = False + for j, B in enumerate(alignments): + if i == j: + continue + if A["qspan"] < B["qspan"] and is_fully_encompassed(A["q_iv"], B["q_iv"]): + encompassed = True + break + if not encompassed: + survivors.append(A) + if not survivors: + survivors = [max(alignments, key=lambda a: a["div"])] + + # 2. Sort by div descending. + survivors.sort(key=lambda a: -a["div"]) + + # 3. Greedy allocation. + claims = [] + claimed_q_all = [] + for X in survivors: + X_q = X["q_iv"] + X_total_mut = X["div"] * X["qspan"] + mutations_in_overlap = 0.0 + for unique_iv, density in claims: + mutations_in_overlap += density * interval_intersection_length(X_q, unique_iv) + X_unique_iv = interval_difference(X_q, claimed_q_all) + X_unique_len = interval_total(X_unique_iv) + X_unique_mut = max(0.0, X_total_mut - mutations_in_overlap) + if X_unique_len > 0: + X_unique_mut = min(X_unique_mut, float(X_unique_len)) + density = X_unique_mut / X_unique_len + else: + X_unique_mut = 0.0 + density = 0.0 + claims.append((X_unique_iv, density)) + claimed_q_all = merge_intervals(claimed_q_all + X_unique_iv) + + total_unique_mut = sum(d * interval_total(iv) for iv, d in claims) + total_unique_len = interval_total(claimed_q_all) + pid = (1.0 - total_unique_mut / total_unique_len) if total_unique_len > 0 else 0.0 + + # eff_qcov, eff_tcov: union of all alignment spans (use ALL alignments, + # not just survivors -- encompassed ones add nothing new anyway). When + # fill_colinear is on, inner gaps within colinear chains are bridged before + # the union so a single element fragmented into HSPs counts as one span. + if fill_colinear: + q_chains, t_chains = chained_intervals( + alignments, max_gap=max_gap, gap_tol=gap_tol) + eff_qcov = interval_total(merge_intervals(q_chains)) / qlen if qlen > 0 else 0.0 + eff_tcov = interval_total(merge_intervals(t_chains)) / tlen if tlen > 0 else 0.0 + else: + all_q_iv = [iv for a in alignments for iv in a["q_iv"]] + all_t_iv = [iv for a in alignments for iv in a["t_iv"]] + eff_qcov = interval_total(merge_intervals(all_q_iv)) / qlen if qlen > 0 else 0.0 + eff_tcov = interval_total(merge_intervals(all_t_iv)) / tlen if tlen > 0 else 0.0 + return pid, eff_qcov, eff_tcov + + +def process_paf(lines, min_pid=0.70, min_qcov=0.70, min_tcov=0.70, + fill_colinear=False, max_gap=5000, gap_tol=0.20, verbose=False): + """Stream PAF lines, group by (qname,tname), pick best target per query, + return [(qname, pass_str, pid, qcov, tcov, best_tname), ...] sorted by qname. + """ + pair_alns = defaultdict(list) + pair_lengths = {} + n_lines = 0 + n_skipped = 0 + for lineno, raw in enumerate(lines, start=1): + if not raw.strip() or raw.startswith("#") or raw.startswith("["): + continue + rec = parse_paf_line(raw, lineno, strict=False) + if rec is None: + n_skipped += 1 + continue + key = (rec["qname"], rec["tname"]) + pair_alns[key].append(rec) + if key not in pair_lengths: + pair_lengths[key] = (rec["qlen"], rec["tlen"]) + n_lines += 1 + if verbose: + print("[classify_ltr_paf_fast] parsed %d alignments across %d (q,t) pairs" + % (n_lines, len(pair_alns)), file=sys.stderr) + if n_skipped: + print("[classify_ltr_paf_fast] skipped %d malformed lines (missing dv:f/de:f or unparseable)" + % n_skipped, file=sys.stderr) + + per_query = defaultdict(list) + for (qname, tname), alns in pair_alns.items(): + qlen, tlen = pair_lengths[(qname, tname)] + pid, qcov, tcov = compute_pair_metrics( + alns, qlen, tlen, + fill_colinear=fill_colinear, max_gap=max_gap, gap_tol=gap_tol) + passes = (pid >= min_pid) and (qcov >= min_qcov) and (tcov >= min_tcov) + per_query[qname].append({ + "tname": tname, + "pid": pid, + "qcov": qcov, + "tcov": tcov, + "passes": passes, + "joint": pid * qcov * tcov, + }) + + results = [] + n_pass = 0 + for qname in sorted(per_query): + candidates = per_query[qname] + candidates.sort(key=lambda c: (-int(c["passes"]), -c["joint"], c["tname"])) + best = candidates[0] + if best["passes"]: + n_pass += 1 + results.append(( + qname, + "pass" if best["passes"] else "fail", + best["pid"], + best["qcov"], + best["tcov"], + best["tname"], + )) + if verbose: + print("[classify_ltr_paf_fast] %d/%d queries pass at pid>=%.3f qcov>=%.3f tcov>=%.3f" + % (n_pass, len(results), min_pid, min_qcov, min_tcov), file=sys.stderr) + return results + + +def format_row(row): + qname, pass_str, pid, qcov, tcov, tname = row + return "%s\t%s\t%.4f\t%.4f\t%.4f\t%s" % (qname, pass_str, pid, qcov, tcov, tname) + + +def main(): + ap = argparse.ArgumentParser( + description="Classify putative LTR-RTs from a minimap2 PAF (CIGAR-free, " + "uses dv:f or de:f and standard PAF columns only).") + ap.add_argument("paf", help="input PAF (use - for stdin)") + ap.add_argument("-o", "--output", default="-", + help="output TSV path (default: stdout)") + ap.add_argument("--min-pid", type=float, default=0.70) + ap.add_argument("--min-qcov", type=float, default=0.70) + ap.add_argument("--min-tcov", type=float, default=0.70) + ap.add_argument("--fill-colinear-gaps", action="store_true", + help="bridge inner gaps within colinear HSP chains " + "(same strand, q-order matches t-order, q-gap ~ t-gap) " + "before computing eff_qcov/eff_tcov. Off by default.") + ap.add_argument("--bridge-max-gap", type=int, default=5000, + help="max inner gap (bp) to bridge on either axis (default: 5000).") + ap.add_argument("--bridge-gap-tol", type=float, default=0.20, + help="max relative mismatch |q_gap - t_gap| / max(q_gap, t_gap) " + "to treat as a synchronized indel (default: 0.20).") + ap.add_argument("--header", action="store_true") + ap.add_argument("-v", "--verbose", action="store_true") + args = ap.parse_args() + + in_fh = sys.stdin if args.paf == "-" else open(args.paf, "r") + out_fh = sys.stdout if args.output == "-" else open(args.output, "w") + if args.verbose: + print("[classify_ltr_paf_fast] reading %s" % args.paf, file=sys.stderr) + try: + results = process_paf( + in_fh, + min_pid=args.min_pid, min_qcov=args.min_qcov, min_tcov=args.min_tcov, + fill_colinear=args.fill_colinear_gaps, + max_gap=args.bridge_max_gap, gap_tol=args.bridge_gap_tol, + verbose=args.verbose, + ) + finally: + if args.paf != "-": + in_fh.close() + try: + if args.header: + out_fh.write("qname\tpass\tpid\teff_qcov\teff_tcov\tbest_tname\n") + for row in results: + out_fh.write(format_row(row) + "\n") + finally: + if args.output != "-": + out_fh.close() + if args.verbose: + print("[classify_ltr_paf_fast] done", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/tesorter2/minimap.py b/tesorter2/minimap.py new file mode 100644 index 0000000..77ad566 --- /dev/null +++ b/tesorter2/minimap.py @@ -0,0 +1,193 @@ +""" +minimap.py — minimap2 wrapper for pass-2 similarity search. + +Pass-2 uses minimap2 with sensitivity-tuned flags chosen to maximize LTR-RT +recall at moderate identity (benchmarked F1 ~0.90 at the 70-70-70 rule). The +PAF is consumed by classify_ltr_paf_fast.process_paf, which groups chains per +(query, target) pair, picks the best target per query, and emits one TSV row +per query with pass/fail status under the user-supplied I-C-L rule. + +PAFRecord / parse_paf_besthit / besthit_per_query are kept for callers that +still want chain-union best-hit semantics, but the live pass-2 pipeline no +longer uses them. +""" + +import logging +import os +import shutil +import subprocess +from collections import OrderedDict, defaultdict + +log = logging.getLogger(__name__) + + +def check_minimap2(bin="minimap2"): + if shutil.which(bin) is None: + raise RuntimeError( + f"{bin!r} not found on PATH. Install minimap2 " + f"(conda: `mamba install -c bioconda minimap2`) and retry." + ) + + +def minimap2_version(bin="minimap2"): + if shutil.which(bin) is None: + log.warning(f"{bin!r} not found on PATH") + return None + r = subprocess.run([bin, "--version"], capture_output=True, text=True, check=False) + v = (r.stdout or r.stderr).strip().splitlines()[0] if (r.stdout or r.stderr) else "unknown" + log.info(f"minimap2 version: {v}") + return v + + +def run_minimap2(query_fa, target_fa, paf_out, ncpu=4, + preset="asm20", extra="", + minimap2_bin="minimap2"): + """ + Run minimap2 and write PAF. + + Sensitivity-tuned flags for LTR-RT pass-2: + -x {preset} base preset (default asm20 = ~20% divergence) + --rmq=no disable repeat-mask query mode + --no-long-join do not extend chains across long gaps + -k 10 -w 10 smaller k-mer / window than asm20 default + -r 500,20000 chain bandwidth bounds + -g 500 stop chain extension at 500 bp gap + -p 0.3 keep secondaries scoring >=30% of primary + -N 100 up to 100 secondaries per query + -m 30 min chaining score + -K 1G large minibatch for throughput + --seed 11 deterministic seeding + --paf-no-hit emit placeholder PAF lines for unmapped queries + -t ncpu + """ + os.makedirs(os.path.dirname(os.path.abspath(paf_out)) or ".", exist_ok=True) + + cmd = ( + f"{minimap2_bin} -x {preset} --rmq=no --no-long-join " + f"-k 10 -w 10 -r 500,20000 -g 500 -p 0.3 -N 100 -m 30 " + f"-t {ncpu} -K 1G --seed 11 --paf-no-hit " + f"{extra} " + f"-o {paf_out} {target_fa} {query_fa}" + ) + log.info(f"minimap2 cmd: {cmd}") + r = subprocess.run(cmd, shell=True, capture_output=True, text=True, check=False) + if r.returncode != 0: + raise RuntimeError( + f"minimap2 failed (exit {r.returncode}): " + f"{(r.stderr or r.stdout)[:2000]}" + ) + return paf_out + + +class PAFRecord: + """One minimap2 PAF line. 0-based half-open coords on query and target.""" + __slots__ = ( + "qseqid", "qlen", "qstart", "qend", "strand", + "sseqid", "tlen", "tstart", "tend", + "matches", "alnlen", "mapq", "score", + ) + + def __init__(self, line): + vals = line.rstrip("\n").split("\t") + self.qseqid = vals[0] + self.qlen = int(vals[1]) + self.qstart = int(vals[2]) + self.qend = int(vals[3]) + self.strand = vals[4] + self.sseqid = vals[5] + self.tlen = int(vals[6]) + self.tstart = int(vals[7]) + self.tend = int(vals[8]) + self.matches = int(vals[9]) + self.alnlen = int(vals[10]) + self.mapq = int(vals[11]) + # AS tag = Smith-Waterman alignment score (when -c is used) + self.score = None + for tag in vals[12:]: + if tag.startswith("AS:i:"): + self.score = int(tag[5:]) + break + if self.score is None: + # fallback: use matches as score proxy + self.score = self.matches + + +def _union_len(intervals): + """Length of the union of [lo, hi) half-open intervals.""" + if not intervals: + return 0 + intervals = sorted(intervals) + total = 0 + cur_lo, cur_hi = intervals[0] + for lo, hi in intervals[1:]: + if lo <= cur_hi: + cur_hi = max(cur_hi, hi) + else: + total += cur_hi - cur_lo + cur_lo, cur_hi = lo, hi + total += cur_hi - cur_lo + return total + + +def parse_paf_besthit(paf_path): + """Parse PAF, union per (query, target) pair on both axes, return one record + per (qseqid, sseqid) pair with merged qcov/tcov/identity; then best-hit-per- + query is selected by caller. + + Returns: list of merged hit dicts with keys: + qseqid, sseqid, qlen, tlen, + qcov (0..1), tcov (0..1), fident (0..1), + alnlen (union on query axis, i.e. unique query bases aligned), + matches, score (max AS across chains for this pair) + """ + if not os.path.exists(paf_path) or os.path.getsize(paf_path) == 0: + return [] + + groups = defaultdict(list) + with open(paf_path) as f: + for line in f: + if not line.strip(): + continue + r = PAFRecord(line) + groups[(r.qseqid, r.sseqid)].append(r) + + merged = [] + for (q, t), recs in groups.items(): + q_intervals = [(r.qstart, r.qend) for r in recs] + t_intervals = [(r.tstart, r.tend) for r in recs] + q_union = _union_len(q_intervals) + t_union = _union_len(t_intervals) + + qlen = recs[0].qlen + tlen = recs[0].tlen + + total_aln = sum(r.alnlen for r in recs) + total_matches = sum(r.matches for r in recs) + fident = (total_matches / total_aln) if total_aln else 0.0 + + best_score = max(r.score for r in recs) + + merged.append({ + "qseqid": q, + "sseqid": t, + "qlen": qlen, + "tlen": tlen, + "qcov": q_union / qlen if qlen else 0.0, + "tcov": t_union / tlen if tlen else 0.0, + "fident": fident, + "alnlen": q_union, + "matches": total_matches, + "score": best_score, + }) + return merged + + +def besthit_per_query(merged): + """Return OrderedDict[qseqid -> best merged record by score].""" + best = OrderedDict() + # Stable sort by score desc for determinism. + merged_sorted = sorted(merged, key=lambda m: -m["score"]) + for m in merged_sorted: + if m["qseqid"] not in best: + best[m["qseqid"]] = m + return best diff --git a/tesorter2/pass2_external.py b/tesorter2/pass2_external.py new file mode 100644 index 0000000..e430631 --- /dev/null +++ b/tesorter2/pass2_external.py @@ -0,0 +1,213 @@ +""" +pass2_external.py — helpers for the --pass2-classified-fasta feature. + +Ported from github.com/cwb14/TEsorter branch `my-new-idea2` (TEsorter/app.py, +head commit b398509). The upstream helpers used Biopython SeqIO and TEsorter's +CommonClassification namedtuple. Here we use pyfastx (already a TEsorter2 +dependency) and emit dicts matching TEsorter2's classifications-dict shape +(`id/order/superfamily/clade/complete/strand/domains/score/secondary`). +""" + +import logging +import os +import re + +import pyfastx + +log = logging.getLogger(__name__) + + +_COORD_HEADER_RE = re.compile( + r'^(?P\S+?:\d+[-\.]+\d+)#(?P[^/]+)/(?P[^/]+)/(?P\S+)$' +) + + +def _format_gff_id(s): + """TEsorter's format_gff_id: strip anything after a '#'. Trivial but kept + as a named helper so the intent reads.""" + return s.split("#", 1)[0] + + +def parse_cls_from_fasta_header(header): + """Parse `>id#Order/Superfamily/Clade` -> (id, order, superfamily, clade). + + Returns None for headers without a '#' or without at least Order/Superfamily. + Missing slots are filled with 'Unknown'/'unknown'. + """ + h = header.strip() + if h.startswith('>'): + h = h[1:] + h = h.split(None, 1)[0] + + if '#' not in h: + return None + raw_id, cls = h.split('#', 1) + raw_id = _format_gff_id(raw_id) + + parts = cls.split('/') + if len(parts) < 2: + return None + + order = parts[0] or 'Unknown' + superfamily = parts[1] or 'unknown' + clade = parts[2] if len(parts) >= 3 and parts[2] else 'unknown' + return raw_id, order, superfamily, clade + + +def extend_hmm_classifications_from_fasta(hmm_cls, fasta_path, db_seq_to_dbs): + """Merge external FASTA classifications into hmm_cls in place. + + Each added entry is shaped like TEsorter2's reconciled classifications + dict so classifier.py:568 still sees the fields it expects. Also extends + db_seq_to_dbs so classify_from_blast accepts hits pointing at these IDs. + """ + if fasta_path is None: + return + + added = 0 + skipped = 0 + # pyfastx is faster than Biopython and already in TEsorter2's deps. + fa = pyfastx.Fasta(fasta_path, build_index=True, full_name=True) + for rec in fa: + parsed = parse_cls_from_fasta_header(rec.name) + if not parsed: + skipped += 1 + continue + sid, order, superfamily, clade = parsed + if sid in hmm_cls: + continue + hmm_cls[sid] = { + "id": sid, + "order": order, + "superfamily": superfamily, + "clade": clade, + "complete": "none", + "strand": "?", + "domains": "none", + "score": 0.0, + "secondary": [], + } + db_seq_to_dbs.setdefault(sid, set()).add("external") + added += 1 + + log.info( + f"extended pass-1 classifications with {added} entries from {fasta_path} " + f"({skipped} headers skipped: not parseable)" + ) + + +def merge_classified_fastas(out_fa, fa_primary, fa_extra=None, clean_nucl=True): + """Write out_fa as the combined pass-2 target FASTA. + + IDs are stripped of any trailing '#...' annotation so they round-trip + cleanly through mmseqs. Primary takes precedence on duplicate IDs. + When clean_nucl=True, non-ATCG characters are stripped before writing. + """ + seen = set() + n = 0 + + def emit_records(path, fout): + nonlocal n + fa = pyfastx.Fasta(path, build_index=True) + for rec in fa: + rid = _format_gff_id(rec.name) + if rid in seen: + continue + seen.add(rid) + seq = str(rec.seq) + if clean_nucl: + seq = "".join(c for c in seq.upper() if c in "ATCG") + fout.write(f">{rid}\n{seq}\n") + n += 1 + + with open(out_fa, "w") as fout: + emit_records(fa_primary, fout) + if fa_extra: + emit_records(fa_extra, fout) + + suffix = " [non-ATCG stripped]" if clean_nucl else "" + log.info(f"pass-2 database FASTA written: {out_fa} ({n} unique IDs){suffix}") + + +def update_classified_fasta_headers(fasta_path, hmm_cls, tmpdir): + """Rewrite a pass2-classified FASTA, upgrading `unknown` slots from hmm_cls. + + Only headers shaped `>chr:start-end#Order/Superfamily/Clade` where the ID + matches an entry in hmm_cls AND at least one of Order/Superfamily/Clade + contains 'unknown' are candidates. Returns the path to the written FASTA, + or None if the input was None. + """ + if fasta_path is None: + return None + + os.makedirs(tmpdir, exist_ok=True) + updated_path = os.path.join(tmpdir, "pass2_classified_updated.fa") + n_updated = 0 + n_total = 0 + + fa = pyfastx.Fasta(fasta_path, build_index=True, full_name=True) + with open(updated_path, "w") as fout: + for rec in fa: + n_total += 1 + header = rec.name.split(None, 1)[0] + seq = str(rec.seq) + + if '#' in header: + raw_id_part, _ = header.split('#', 1) + rid = _format_gff_id(raw_id_part) + + m = _COORD_HEADER_RE.match(header) + if m and rid in hmm_cls: + old_order = m.group('order') + old_sfam = m.group('sfam') + old_clade = m.group('clade') + + if 'unknown' in (old_order.lower(), old_sfam.lower(), old_clade.lower()): + cls = hmm_cls[rid] + new_cls = "{}/{}/{}".format( + cls["order"], cls["superfamily"], cls["clade"] + ) + header = f"{raw_id_part}#{new_cls}" + n_updated += 1 + + fout.write(f">{header}\n{seq}\n") + + log.info( + f"updated {n_updated}/{n_total} headers in pass2-classified-fasta " + f"using pass-1 classifications" + ) + return updated_path + + +def clean_fasta_atcg(path): + """In-place ATCG-only cleaner. Kept for completeness / upstream parity. + + Not actively called by TEsorter2's pass-2 because the mmseqs wrapper + already cleans the query and merge_classified_fastas cleans the DB. + """ + tmp = path + ".atcg_clean.tmp" + with open(path) as fin, open(tmp, "w") as fout: + buf = [] + header = [None] + + def flush(): + if header[0] is None: + return + seq = "".join(buf).upper() + seq = "".join(c for c in seq if c in "ATCG") + fout.write(header[0] + "\n") + fout.write(seq + "\n") + + for line in fin: + line = line.rstrip("\n") + if not line: + continue + if line.startswith(">"): + flush() + header[0] = line + buf.clear() + else: + buf.append(line) + flush() + os.replace(tmp, path) + log.info(f"non-ATCG characters removed from {path}") diff --git a/tesorter2/pipeline.py b/tesorter2/pipeline.py index b5708df..4d105d1 100644 --- a/tesorter2/pipeline.py +++ b/tesorter2/pipeline.py @@ -2,7 +2,8 @@ Main pipeline for TE classification. Orchestrates: FASTA ingestion -> alphabet detection -> optional translation --> HMM search -> classification -> BLAST pass-2 -> SQLite + TSV output. +-> HMM search -> classification -> pass-2 similarity search -> SQLite + TSV +output. """ import argparse @@ -25,6 +26,7 @@ store_classifications, reconcile_classifications, DB_CONFIGS) from .blast_pass2 import blast_pass2 +from .minimap import minimap2_version from . import bath_search @@ -210,6 +212,56 @@ def parse_args(): "per-database classifications and their summed normalized " "scores in descending order of evidence strength.", ) + + # pass-2 options + parser.add_argument( + "-dp2", "--disable-pass2", + action="store_true", default=False, + help="Skip pass-2 similarity search (HMM-only classification)", + ) + parser.add_argument( + "-rule", "--pass2-rule", + default="80-80-80", type=str, metavar="I-C-L", + help="Pass-2 threshold as identity-coverage-length. For the blast " + "aligner: pident, qcovs, and alignment-length filters (80-80-80 " + "matches TEsorter2 master). For minimap2: I drives " + "classify_ltr_paf_fast --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 [default: %(default)s]", + ) + parser.add_argument( + "--pass2-classified-fasta", + default=None, type=str, metavar="FASTA", + help="Optional FASTA of previously-classified elements to augment " + "the pass-2 target database. Headers must be like " + ">id#Order/Superfamily/Clade", + ) + parser.add_argument( + "--minimap2-extra", + default="", type=str, metavar="STR", + help="Extra flags passed through to minimap2 (advanced) " + "[default: empty]", + ) + parser.add_argument( + "--pass2-aligner", + choices=["blast", "minimap2"], default="blast", + help="Aligner for the pass-2 similarity search. 'blast' (default) " + "reproduces TEsorter2 master's blastn pass-2 (qcovs + " + "alignment-length filter, clade=unknown); 'minimap2' uses the " + "PAF qcov+tcov path and inherits the best target's full " + "classification. Both share the same -rule and the " + "--pass2-classified-fasta external-pool merge.", + ) + parser.add_argument( + "--blast-task", + choices=["megablast", "dc-megablast"], default="megablast", + help="blastn -task for the 'blast' pass-2 aligner. megablast " + "(default) is fastest and tuned for near-identical matches; " + "dc-megablast uses discontiguous seeds — slower but more " + "sensitive to diverged/cross-species matches. Ignored when " + "--pass2-aligner=minimap2 [default: %(default)s]", + ) + parser.add_argument( "--no-tesorter-outputs", action="store_true", @@ -574,16 +626,37 @@ def main(): log.info(f" Reconciled across {len(per_db_results)} databases: " f"{len(reconciled)} sequences") - # --- BLAST pass-2 --- + # --- pass-2 similarity search --- all_results = list(reconciled) - if not args.pass_1_only and all_classifications: - log.info("--- BLAST pass-2 ---") + if (not args.pass_1_only and not args.disable_pass2 + and all_classifications): + try: + p2_id, p2_cov, p2_len = args.pass2_rule.split("-") + p2_id = float(p2_id) + p2_cov = float(p2_cov) + p2_len = float(p2_len) + except ValueError: + raise SystemExit( + f"--pass2-rule must be I-C-L (three numbers separated by '-'), " + f"got {args.pass2_rule!r}" + ) + + log.info(f"--- pass-2 ({args.pass2_aligner}) ---") + if args.pass2_aligner == "minimap2": + minimap2_version() blast_cls = blast_pass2( args.sequence, conn, hmm_classifications=all_classifications, seq_type="nucl", n_processors=args.processors, + min_identity=p2_id, + min_coverage=p2_cov, + min_length=p2_len, outdir=outdir, + pass2_classified_fasta=args.pass2_classified_fasta, + minimap2_extra=args.minimap2_extra, + aligner=args.pass2_aligner, + blast_task=args.blast_task, ) if blast_cls: diff --git a/tesorter2/tesorter_compat.py b/tesorter2/tesorter_compat.py index cb37ba5..88da424 100644 --- a/tesorter2/tesorter_compat.py +++ b/tesorter2/tesorter_compat.py @@ -71,7 +71,29 @@ def parse_args(): help="Minimum normalized score [default: 0.1]") parser.add_argument("-dp2", "--disable-pass2", action="store_true", default=False, - help="Do not run BLAST pass-2 classification") + help="Do not run pass-2 classification") + parser.add_argument("-rule", "--pass2-rule", type=str, default="80-80-80", + metavar="I-C-L", + help="Pass-2 threshold identity-coverage-length. " + "blast: pident, qcovs, and alignment-length " + "filters. minimap2: I drives " + "classify_ltr_paf_fast --min-pid; C drives both " + "--min-qcov and --min-tcov; L unused " + "[default: 80-80-80]") + parser.add_argument("--pass2-aligner", choices=["blast", "minimap2"], + default="blast", + help="Aligner for the pass-2 similarity search. " + "blast (default) reproduces TEsorter2 master's " + "blastn pass-2; minimap2 uses the PAF qcov+tcov " + "path [default: blast]") + parser.add_argument("--pass2-classified-fasta", type=str, default=None, + metavar="FASTA", + help="Optional FASTA of previously-classified elements " + "to augment pass-2 target DB. Headers must be " + "like >id#Order/Superfamily/Clade") + parser.add_argument("--minimap2-extra", type=str, default="", + metavar="STR", + help="Extra flags passed through to minimap2") parser.add_argument("-nolib", "--no-library", action="store_true", default=False, help="Do not generate RepeatMasker library file") @@ -234,15 +256,30 @@ def main(): export_classification_tsv(results, cls_out) log.info(f"Classification: {len(results)} sequences -> {cls_out}") - # BLAST pass-2 + # pass-2 similarity search if not args.disable_pass2 and args.seq_type == "nucl": + try: + p2_id, p2_cov, p2_len = args.pass2_rule.split("-") + p2_id = float(p2_id) + p2_cov = float(p2_cov) + p2_len = float(p2_len) + except ValueError: + log.error(f"--pass2-rule must be I-C-L, got {args.pass2_rule!r}") + sys.exit(1) + hmm_cls = {r["id"]: r for r in results} blast_cls = blast_pass2( args.sequence, conn, hmm_classifications=hmm_cls, seq_type="nucl", n_processors=args.processors, + min_identity=p2_id, + min_coverage=p2_cov, + min_length=p2_len, outdir=args.tmp_dir or os.path.dirname(prefix) or ".", + pass2_classified_fasta=args.pass2_classified_fasta, + minimap2_extra=args.minimap2_extra, + aligner=args.pass2_aligner, ) if blast_cls: @@ -250,7 +287,7 @@ def main(): mode=run_mode) all_results = results + blast_cls export_classification_tsv(all_results, cls_out) - log.info(f"BLAST pass-2: {len(blast_cls)} additional -> {cls_out}") + log.info(f"minimap2 pass-2: {len(blast_cls)} additional -> {cls_out}") # Generate TEsorter-format output files if config and results: