45:["$","title",null,{"children":"Closest phages and genome checks | WabushOnce | Observable"}] 4d:T3ee57,{ // ========================================================== // Closest PhagesDB phages for a folder of genomes — script generator for Windows (WSL) and Mac (Observable) // Enter a folder with one subfolder per phage (each with a .fasta), and the PhagesDB downloads; then download // the script and run it in WSL. Writes closest_phages.html (linked, sortable) and closest_phages.tsv. // ========================================================== // host genera offered for the GenBank search (true = ticked by default) const GENERA = [["Mycobacterium", true], ["Gordonia", true], ["Microbacterium", true], ["Arthrobacter", true], ["Staphylococcus", true], ["Streptomyces", false], ["Rhodococcus", false], ["Corynebacterium", false], ["Enterococcus", false], ["Streptococcus", false], ["Lactococcus", false], ["Bacillus", false], ["Listeria", false], ["Pseudomonas", false], ["Escherichia", false], ["Klebsiella", false], ["Salmonella", false], ["Acinetobacter", false]]; const wrap = document.createElement("div"); wrap.style.maxWidth = "980px"; wrap.style.fontFamily = "system-ui, -apple-system, Segoe UI, sans-serif"; wrap.innerHTML = `

Closest phages and genome checks — script generator (Windows with WSL, or Mac)

Finds the closest PhagesDB phage for every genome in a folder (one subfolder per phage, each with a .fasta file), using minimap2 against the PhagesDB genome collection. Writes closest_phages.html (sortable table with PhagesDB and GenBank links), closest_phages.tsv (one row per phage) and closest_phages_top_hits.tsv (the best hits for each phage). Optional genome checks: host DNA %, CheckV completeness/contamination and a similarity matrix of all your genomes (similarity_matrix.tsv + heatmap). Works on the Nanopore pipeline's output folder (then each phage links to its coverage.html and shows the pipeline's warnings). Report: status of each genome against its closest phage (complete / fragmented / partial), ICTV-style similarity, same-phage duplicates, filters by host and cluster, and a filterable similarity heatmap.

0. This computer

1. Genomes

Folder with one subfolder per phage

FASTA file in each subfolder — optional

Automatic: <phage>.fasta (the whole assembly), else <phage>.longest_contig.fasta, else the largest FASTA (Flye/Autocycler drafts are skipped). Or a name such as {phage}.fasta or *_final.fasta. FASTA files lying directly in the folder are also searched (named by file).

Number of hits listed per phage

Search every contig of at least (bp)

All contigs this long are searched together, so a genome in pieces still finds its relative; the report shows how much of the relative they cover (complete / fragmented / partial). The longest contig is used when none is this long.

Isolation hosts — optional (for the Host filter and host checks)

Paste phage and host, one per line

or a CSV / TSV file (e.g. a PhageTracker export)

Columns named like phage / name and host / isolation host are found automatically; without a header the first two columns are used. The report's Host filter then uses the isolation host, and phages whose closest relative infects another genus are marked.

2. PhagesDB downloads

PhagesDB genomes FASTA

PhagesDB data table (tab-delimited) — optional

Adds subcluster, host strain, temperate, morphotype, ends and GenBank accession. A name without an extension also finds .txt, .tsv or .tab. The index is built once per downloaded file in a phagesdb_index folder next to the FASTA and shared with the Nanopore phage pipeline.

3. GenBank phages (NCBI) — an extra search, BLAST-style columns

PhagesDB is always searched, with all of its phages (every host: Mycobacterium, Gordonia, Microbacterium, Arthrobacter, Streptomyces …). The choice below only adds a second search of GenBank, for phages that are not in PhagesDB (e.g. Staphylococcus phages, or phages from labs outside SEA-PHAGES).

Also search GenBank (in addition to PhagesDB)

Host genera for the GenBank search

${GENERA.map(([g, on]) => ` ${g}`).join("")}

Tick the host genera of your phages. Only limits the GenBank download; it doesn't affect the PhagesDB search. Phages are kept when their host is any species of these genera or their name starts with one (e.g. "Staphylococcus phage CF9"). The first run downloads all complete tailed-phage genomes once (roughly 0.5–1 GB, several minutes) and keeps only the matching ones.

NCBI genomes FASTA you downloaded yourself — optional

Download the NCBI genomes again now (this also happens automatically when they are more than 90 days old)

Downloaded once (needs internet in WSL) with NCBI's Datasets tool, which the script fetches itself if it isn't installed, into an NCBI_phages folder next to the PhagesDB FASTA, then reused. GenBank matches compete with PhagesDB for the report's "Closest phage" (identity and query cover computed with minimap2, so close to BLASTn's numbers, without E-values); the top GenBank hits are in closest_genbank_top_hits.tsv.

4. Genome checks: host DNA, completeness, similarity

Host genome FASTA files — optional (up to 5)

Each genome gets a Host DNA %: the share that matches a host genome (≥85% identity). Over 50% means the contig is mostly host DNA, not phage. Use the same host FASTAs as in the Nanopore pipeline.

Run CheckV (completeness, contamination, quality) CheckV conda environment Create it if it is missing (conda/mamba, bioconda)

CheckV database folder (WSL) Download the database if it is missing (once; a few GB)

A WSL folder (e.g. ~/databases/checkv) is faster than a Windows one.

Similarity matrix of all genomes

VIRIDIC-style intergenomic similarity between every pair, with species (≥95%) and genus (≥70%) groups (ICTV thresholds) and a heatmap; shows duplicates and which phages belong together.

5. Output and environment

Report folder — optional

Conda environment with minimap2

The Medaka environment from the pipeline has minimap2.

Conda folder — optional

Threads

Open the HTML report when finished

6. Run it

The script carries the search program inside it, so download it and run the file rather than pasting it:

1. Download closest_phages.sh

Not sure the right genomes will be used? Run the check command first: it lists every phage, the FASTA chosen, its contigs and isolation host, without searching.

Paste that line into WSL. It runs the newest closest_phages*.sh in your Downloads folder (also a renamed copy such as closest_phages (1).sh).

`; const $ = (id) => wrap.querySelector(`#${id}`); // v356: status and cluster from each genome's own contigs, completeness evidence, contamination and Ligase-pair notes, clusters table; v355: contig-cluster filters and hatching, colour legend; v354: more-reads check, host DNA and contig-cluster figures, phage names on the read charts; v353: status reasons and grouped table columns; v352: Mac or Windows (the pipeline generator remembers the choice in this browser) try { const pp = localStorage.getItem("phagePlatform"); if (pp && wrap.querySelector("#cpPlatform option[value='" + pp + "']")) wrap.querySelector("#cpPlatform").value = pp; } catch (_) {} const isMac = () => $("cpPlatform").value !== "wsl"; function macPath(value) { let p = String(value ?? "").trim(); if ((p.startsWith('"') && p.endsWith('"')) || (p.startsWith("'") && p.endsWith("'"))) p = p.slice(1, -1).trim(); if (/^file:\/\//i.test(p)) { try { p = decodeURIComponent(new URL(p).pathname); } catch (_) {} } p = p.replace(/\\ /g, " "); return p.length > 1 ? p.replace(/\/+$/, "") : p; } const toPath = (v) => (isMac() ? macPath(v) : toWslPath(v)); function applyPlatformText() { const mac = isMac(), pairs = mac ? [[/needs internet in WSL/g, "needs internet"], [/CheckV database folder $WSL$/g, "CheckV database folder"], [/A WSL folder \(e\.g\./g, "A local folder (e.g."], [/ is faster than a Windows one\./g, "."], [/Paste that line into WSL\./g, "Paste that line into Terminal."]] : []; const tw = document.createTreeWalker(wrap, NodeFilter.SHOW_TEXT); for (let n = tw.nextNode(); n; n = tw.nextNode()) { if (n.parentElement && n.parentElement.closest("pre, script, style, select")) continue; if (n.__osOrig == null) { if (!/WSL|Windows one/.test(n.nodeValue)) continue; n.__osOrig = n.nodeValue; } let t = n.__osOrig; for (const [x, y] of pairs) t = t.replace(x, y); if (n.nodeValue !== t) n.nodeValue = t; } wrap.querySelectorAll("input[placeholder]").forEach((el) => { if (el.dataset.osPh == null) { if (!/[A-Za-z]:\\/.test(el.placeholder)) return; el.dataset.osPh = el.placeholder; } el.placeholder = mac ? el.dataset.osPh.replace(/^[A-Za-z]:\\Users\\[^\\]+\\/, "~/").replace(/\\/g, "/") : el.dataset.osPh; }); $("cpPlatNote").innerHTML = mac ? "Paste the run command into Terminal. minimap2 comes from a conda environment (the pipeline's setup script installs it in reads-assembly-core and medaka-polish); " + "the run command keeps the Mac awake while it works. Folder paths look like ~/Nanopore/Results (drag a folder into Terminal to see its path)." : "The script runs in WSL (Ubuntu). Windows paths such as C:\\Users\\you\\… work: they are read as /mnt/c/Users/you/…."; } let downloadUrl = null; const shQuote = (value) => "'" + String(value).replace(/'/g, "'\"'\"'") + "'"; function toWslPath(value) { let p = String(value ?? "").trim(); if (!p) return ""; if ((p.startsWith('"') && p.endsWith('"')) || (p.startsWith("'") && p.endsWith("'"))) p = p.slice(1, -1).trim(); if (/^file:\/\//i.test(p)) { try { const u = new URL(p); p = decodeURIComponent(u.pathname); if (/^\/[A-Za-z]:\//.test(p)) p = p.slice(1); } catch (_) {} } p = p.replace(/^\\\\\?\\/, ""); const m = p.match(/^([A-Za-z]):[\\/]?(.*)$/); if (m) { const rest = m[2].replace(/\\/g, "/").replace(/^\/+/, "").replace(/\/+$/, ""); return rest ? `/mnt/${m[1].toLowerCase()}/${rest}` : `/mnt/${m[1].toLowerCase()}`; } if (/^\/mnt\/[A-Za-z](\/|$)/.test(p)) p = p.replace(/^\/mnt\/([A-Za-z])/, (_, d) => `/mnt/${d.toLowerCase()}`); return p.length > 1 ? p.replace(/\/+$/, "") : p; } const CLOSEST_PY = `#!/usr/bin/env python3 """ Closest PhagesDB phages (and closest GenBank phages) for a folder of phage genomes. For every subfolder of --input (one per phage) the genome FASTA is found, and all genomes are compared at once with minimap2 (asm20) to the PhagesDB genome collection (Actinobacteriophages-All.fasta). Writes closest_phages.html (sortable, linked table), closest_phages.tsv (one row per phage) and closest_phages_top_hits.tsv (the best hits for every phage). The PhagesDB index is shared with the Nanopore phage pipeline (folder "phagesdb_index" next to the FASTA), so it is built only once per downloaded PhagesDB file. Optionally, complete phage genomes are downloaded once from NCBI (NCBI Datasets: by host genus, or all Caudoviricetes) and searched the same way, giving BLAST-style columns for the closest GenBank phage (description, accession, per. ident, query cover, aligned length, subject length). Genome checks (optional): the share of each genome that matches the host genome(s) (host DNA), CheckV completeness, contamination and quality, and a VIRIDIC-style intergenomic similarity matrix of all the genomes with species (>= 95%) and genus (>= 70%) groups (ICTV thresholds), shown as a heatmap. v2 (2026-11-02): every contig >= --min-contig (5 kb) of a genome is searched (joined with N spacers), not only the longest; the report shows how much of the closest phage your contigs cover (complete / fragmented / partial / no relative), near-equal matches ("ties"), the ICTV-style similarity to the closest phage (same species / same genus / different genus), phages in the run that are the same phage (>= 98% similar), isolation hosts (--isolation-hosts), links to the pipeline's coverage.html and its warnings. --check lists the genomes and FASTA files without searching. """ import argparse import gzip import html import json import os import re import subprocess import sys import time import urllib.parse FASTA_EXT = (".fasta", ".fa", ".fna", ".fas", ".fasta.gz", ".fa.gz", ".fna.gz") TOOL_VERSION = "closest_phages v3.1 (2026-11-06)" SPACER = "N" * 100 MAX_CONTIG = 1_000_000 # longer contigs are bacterial (host) DNA, not phage: not searched # between the contigs of one genome in the query FASTA def open_text(path): return gzip.open(path, "rt") if path.endswith(".gz") else open(path, errors="replace") def read_fasta(path): recs, name, desc, chunks = [], None, "", [] with open_text(path) as fh: for line in fh: line = line.rstrip("\\r\\n") if line.startswith(">"): if name is not None: recs.append((name, desc, "".join(chunks))) parts = line[1:].strip().split(None, 1) name = parts[0] if parts else "contig" desc = parts[1] if len(parts) > 1 else "" chunks = [] elif line.strip(): chunks.append(line.strip()) if name is not None: recs.append((name, desc, "".join(chunks))) return recs def write_fasta(fh, name, seq, width=80): fh.write(">" + name + "\\n") for i in range(0, len(seq), width): fh.write(seq[i:i + width] + "\\n") def merge(iv): out = [] for s, e in sorted(iv): if out and s <= out[-1][1]: out[-1][1] = max(out[-1][1], e) else: out.append([s, e]) return out # ---------------------------------------------------------------- PhagesDB header + data table def pdb_header(header): h = header.strip() out = {"host": "", "phage": "", "length": "", "overhang": "", "ends": "", "cluster": ""} m = re.match(r"(.+?)\\s+phage\\s+([^\\s,;]+)", h, re.I) if m: out["host"], out["phage"] = m.group(1).strip(), m.group(2).strip() else: out["phage"] = h.split(None, 1)[0] if h else "" m = re.search(r"([\\d,]+)\\s*bp", h) if m: out["length"] = m.group(1).replace(",", "") m = re.search(r"(\\d+)[- ]base\\s+([35])['’′]?\\s*overhang\\s*\$([A-Za-z]+)\$", h, re.I) if m: out["overhang"] = m.group(3).upper() out["ends"] = f"{m.group(1)}-base {m.group(2)}' overhang ({m.group(3).upper()})" elif re.search(r"circularly permuted", h, re.I): out["ends"] = "circularly permuted" elif re.search(r"terminal repeat", h, re.I): out["ends"] = "direct terminal repeats" m = re.search(r"Cluster\\s+([A-Za-z0-9]+)", h) if m: out["cluster"] = m.group(1) elif re.search(r"\\bSingleton\\b", h, re.I): out["cluster"] = "Singleton" return out TABLE_COLS = ["in_table", "subcluster", "temperate", "gc_pct", "end_type", "overhang_len", "term_rpt_len", "orfs", "trnas", "morphotype", "host_strain", "accession", "genome_length_table"] def none_(v): v = str(v or "").strip() return "" if v.lower() in ("none", "null", "nan", "n/a") else v def read_table(path): rows = {} if not path or not os.path.isfile(path): return rows raw = open(path, "rb").read() text = "" for enc in ("utf-8-sig", "cp1252", "latin-1"): try: text = raw.decode(enc) break except UnicodeDecodeError: continue lines = text.splitlines() if not lines: return rows cols = [c.strip() for c in lines[0].split("\\t")] for line in lines[1:]: f = line.split("\\t") if not f or not f[0].strip(): continue d = {c: (f[i].strip() if i < len(f) else "") for i, c in enumerate(cols)} rows[d.get("Phage Name", f[0]).strip().lower()] = d return rows def prepare_phagesdb(fasta, table, threads): """IDs FASTA, metadata table and minimap2 index, in phagesdb_index next to the FASTA (same files as the Nanopore phage pipeline uses, so either one builds them once).""" pdb_dir = os.path.join(os.path.dirname(os.path.abspath(fasta)), "phagesdb_index") os.makedirs(pdb_dir, exist_ok=True) st = os.stat(fasta) stem = os.path.join(pdb_dir, f"{os.path.basename(fasta)}.{st.st_size}_{int(st.st_mtime)}.v9") tsig = "none" if table: ts = os.stat(table) tsig = f"{ts.st_size}_{int(ts.st_mtime)}" ids, meta, index = stem + ".ids.fasta", f"{stem}.table_{tsig}.meta.tsv", stem + ".asm20.mmi" if not (os.path.isfile(ids) and os.path.isfile(meta) and os.path.getsize(meta) > 0): print("Reading the PhagesDB genomes" + (" and data table" if table else "") + " (once per downloaded file)…", flush=True) tab = read_table(table) used, n, matched = set(), 0, 0 with open_text(fasta) as fin, open(ids + ".part", "w") as fo, open(meta + ".part", "w") as fm: fm.write("id\\tphage\\thost\\tlength_bp\\tends\\toverhang\\tcluster\\theader\\t" + "\\t".join(TABLE_COLS) + "\\n") for line in fin: if line.startswith(">"): header = line[1:].strip() f = pdb_header(header) base = re.sub(r"[^A-Za-z0-9._-]+", "_", f["phage"] or "phage").strip("_") or "phage" pid, k = base, 1 while pid in used: k += 1 pid = f"{base}_{k}" used.add(pid) n += 1 fo.write(f">{pid}\\n") extra = {c: "" for c in TABLE_COLS} row = tab.get(f["phage"].lower()) if row: matched += 1 sub, clu = none_(row.get("Subcluster")), none_(row.get("Cluster")) f["cluster"] = sub or clu or f["cluster"] end_type = none_(row.get("End Type")) oh_len, oh_seq = none_(row.get("Overhang Length (bp)")), none_(row.get("Overhang Sequence")).upper() trl = none_(row.get("Term Rpt Length")) if oh_seq and re.fullmatch(r"[ACGT]+", oh_seq): f["overhang"] = oh_seq if end_type: bits = [end_type] + ([f"{oh_len}-base"] if oh_len else []) + \\ ([f"({f['overhang']})"] if f["overhang"] else []) + ([f"terminal repeat {trl} bp"] if trl else []) f["ends"] = " ".join(bits) if none_(row.get("Host")): f["host"] = none_(row.get("Host")) extra.update({"in_table": "yes", "subcluster": sub or clu, "temperate": none_(row.get("Temperate?")), "gc_pct": none_(row.get("GC%")), "end_type": end_type, "overhang_len": oh_len, "term_rpt_len": trl, "orfs": none_(row.get("# ORFs")), "trnas": none_(row.get("# tRNAs")), "morphotype": none_(row.get("Morphotype")), "host_strain": none_(row.get("Host")), "accession": none_(row.get("Accession #")), "genome_length_table": none_(row.get("Genome Length(bp)"))}) clean = header.replace("\\t", " ") fm.write(f"{pid}\\t{f['phage']}\\t{f['host']}\\t{f['length']}\\t{f['ends']}\\t{f['overhang']}\\t{f['cluster']}" f"\\t{clean}\\t" + "\\t".join(extra[c].replace("\\t", " ") for c in TABLE_COLS) + "\\n") elif line.strip(): fo.write(line.strip() + "\\n") if n == 0: sys.exit("No sequences in the PhagesDB FASTA.") os.replace(ids + ".part", ids) os.replace(meta + ".part", meta) print(f" {n} PhagesDB genomes" + (f"; {matched} matched to the data table ({len(tab)} rows)" if tab else "")) if not os.path.isfile(index): print("Indexing the PhagesDB genomes for minimap2 (once per downloaded file; a few minutes)…", flush=True) r = subprocess.run(["minimap2", "-x", "asm20", "-t", str(threads), "-d", index + ".part", ids], capture_output=True, text=True) if r.returncode != 0: sys.exit("minimap2 could not index the PhagesDB genomes:\\n" + r.stderr[-800:]) os.replace(index + ".part", index) m = {} with open(meta) as fh: cols = next(fh).rstrip("\\n").split("\\t") for line in fh: f = line.rstrip("\\n").split("\\t") m[f[0]] = dict(zip(cols, f)) return index, m # ---------------------------------------------------------------- NCBI phage genomes def _datasets_url(): """the NCBI Datasets tool for this computer (Linux / WSL, or macOS)""" import platform base = "https://ftp.ncbi.nlm.nih.gov/pub/datasets/command-line/v2/" if platform.system() == "Darwin": return base + "mac/datasets" return base + ("linux-arm64/datasets" if platform.machine().lower() in ("aarch64", "arm64") else "linux-amd64/datasets") DATASETS_URL = _datasets_url() def datasets_tool(ncbi_dir): """NCBI Datasets command-line tool: on PATH, or downloaded once next to the NCBI genomes.""" import shutil import urllib.request found = shutil.which("datasets") if found: return found local = os.path.join(ncbi_dir, "datasets") if os.path.isfile(local) and os.access(local, os.X_OK): return local print(f"Downloading the NCBI Datasets tool → {local}", flush=True) urllib.request.urlretrieve(DATASETS_URL, local + ".part") os.chmod(local + ".part", 0o755) os.replace(local + ".part", local) return local def ncbi_download(ncbi_dir, mode, hosts, stem): """Complete phage genomes from NCBI into STEM.fasta + STEM.meta.tsv (accession, description, host, length). Host-genus mode: NCBI's --host filter only matches the exact host taxon (a phage of "Staphylococcus epidermidis" is not returned for "Staphylococcus") and skips records without a host, so all complete Caudoviricetes genomes are downloaded once and kept when the host lineage includes one of the genera, the host name starts with one, or the title starts with it ("Staphylococcus phage CF9", "... virus ..."). The per-host downloads are added too, in case a matching phage is not classified as Caudoviricetes.""" import zipfile tool = datasets_tool(ncbi_dir) genera = {h.lower() for h in hosts} if mode == "all": jobs = [["Caudoviricetes", None, False]] else: jobs = [["Caudoviricetes", None, True]] + [["10239", h, False] for h in hosts] # 10239 = Viruses def wanted(desc, host_name, lineage): if any(x in genera for x in lineage): return True if host_name and host_name.split()[0].lower() in genera: return True m = re.match(r"\\s*(?:uncultured\\s+)?([A-Za-z]+)\\s+(?:phage|virus|bacteriophage)\\b", desc or "", re.I) return bool(m and m.group(1).lower() in genera) seen, n = set(), 0 with open(stem + ".fasta.part", "w") as fo, open(stem + ".meta.tsv.part", "w") as fm: fm.write("accession\\tdescription\\thost\\tlength_bp\\n") for taxon, host, filt in jobs: label = f"viruses with host {host}" if host else (f"{taxon} (kept: phages of {', '.join(hosts)})" if filt else taxon) print(f" NCBI Datasets: complete genomes of {label} …" + (" (large, once: roughly 0.5–1 GB; the matching phages are kept)" if taxon == "Caudoviricetes" else ""), flush=True) # Big downloads go to a Linux folder (writing them to /mnt/c can fail); a Caudoviricetes package you # downloaded yourself can be reused by saving it as NCBI_phages/Caudoviricetes_complete.zip. tmp_dir = os.path.join(os.path.expanduser("~"), ".cache", "closest_phages") os.makedirs(tmp_dir, exist_ok=True) reuse = os.path.join(ncbi_dir, "Caudoviricetes_complete.zip") if taxon == "Caudoviricetes" else "" keep_zip = bool(reuse) and os.path.isfile(reuse) and (time.time() - os.path.getmtime(reuse)) / 86400 <= NCBI_MAX_AGE_DAYS zpath = reuse if keep_zip else os.path.join(tmp_dir, f"ncbi_{os.getpid()}_{re.sub(r'[^A-Za-z0-9]+', '_', label)[:40]}.zip") if keep_zip: print(f" using {reuse} (delete it to download a newer package)", flush=True) else: cmd = [tool, "download", "virus", "genome", "taxon", taxon, "--complete-only", "--include", "genome", "--filename", zpath] if host: cmd += ["--host", host] r = subprocess.run(cmd, capture_output=True, text=True) if r.returncode != 0 or not os.path.isfile(zpath): msg = (r.stderr or r.stdout or "").strip()[-800:] if os.path.exists(zpath): os.remove(zpath) if taxon == "Caudoviricetes": # without it the GenBank set would be badly incomplete, so stop the GenBank search here raise RuntimeError("NCBI Datasets could not download the Caudoviricetes genomes:\\n" + msg + "\\nTip: download it yourself with\\n datasets download virus genome taxon Caudoviricetes " f"--complete-only --include genome --filename '{reuse or 'Caudoviricetes_complete.zip'}'\\n" "and run the script again.") print(f"⚠️ NCBI Datasets failed for {label} (skipped):\\n{msg}", flush=True) continue hostmap, lineages = {}, {} with zipfile.ZipFile(zpath) as z: names = z.namelist() rep_name = next((x for x in names if x.endswith("data_report.jsonl")), None) if rep_name: with z.open(rep_name) as rh: for raw in rh: try: d = json.loads(raw) except ValueError: continue acc = d.get("accession", "") hd = d.get("host") or {} h = hd.get("organismName") or hd.get("organism_name") or hd.get("name") or "" if acc: hostmap[acc] = h lineages[acc] = {str(x.get("name", "")).lower() for x in (hd.get("lineage") or []) if isinstance(x, dict)} fna = next((x for x in names if x.endswith(".fna")), None) added = 0 if fna: acc, desc, chunks, keep = None, "", [], False def flush(): nonlocal n, added if acc and keep: seq = "".join(chunks) write_fasta(fo, acc, seq) fm.write(f"{acc}\\t{desc.replace(chr(9), ' ')}\\t{hostmap.get(acc, host or '')}\\t{len(seq)}\\n") n += 1 added += 1 with z.open(fna) as fh: for raw in fh: line = raw.decode("utf-8", "replace").rstrip("\\r\\n") if line.startswith(">"): flush() parts = line[1:].split(None, 1) acc = parts[0] if parts else "" desc = parts[1] if len(parts) > 1 else "" keep = acc not in seen and (not filt or wanted(desc, hostmap.get(acc, ""), lineages.get(acc, set()))) if keep: seen.add(acc) chunks = [] elif keep: chunks.append(line.strip()) flush() if not keep_zip: os.remove(zpath) print(f" {added:,} genomes", flush=True) if n == 0: raise RuntimeError("NCBI returned no genomes; check the host names") os.replace(stem + ".fasta.part", stem + ".fasta") os.replace(stem + ".meta.tsv.part", stem + ".meta.tsv") with open(stem + ".info.txt", "w") as fh: fh.write(f"downloaded\\t{time.strftime('%Y-%m-%d %H:%M')}\\nmode\\t{mode}\\nhosts\\t{' '.join(hosts)}\\ngenomes\\t{n}\\n") NCBI_MAX_AGE_DAYS = 90 def ncbi_age_days(stem): try: with open(stem + ".info.txt") as fh: for line in fh: k, _, v = line.rstrip("\\n").partition("\\t") if k == "downloaded": return (time.time() - time.mktime(time.strptime(v.strip(), "%Y-%m-%d %H:%M"))) / 86400 except (OSError, ValueError): pass return None def prepare_ncbi(ncbi_dir, mode, hosts, local_fasta, refresh, threads): """Index of NCBI phage genomes (downloaded once, or a FASTA you downloaded yourself) and its metadata.""" os.makedirs(ncbi_dir, exist_ok=True) if local_fasta: st = os.stat(local_fasta) stem = os.path.join(ncbi_dir, f"{os.path.basename(local_fasta)}.{st.st_size}_{int(st.st_mtime)}") if not os.path.isfile(stem + ".meta.tsv"): with open_text(local_fasta) as fin, open(stem + ".fasta", "w") as fo, open(stem + ".meta.tsv", "w") as fm: fm.write("accession\\tdescription\\thost\\tlength_bp\\n") acc, desc, chunks = None, "", [] for line in list(fin) + [">"]: line = line.rstrip("\\r\\n") if line.startswith(">"): if acc: seq = "".join(chunks) write_fasta(fo, acc, seq) fm.write(f"{acc}\\t{desc}\\t\\t{len(seq)}\\n") parts = line[1:].split(None, 1) acc = parts[0] if parts else None desc = (parts[1] if len(parts) > 1 else "").replace("\\t", " ") chunks = [] else: chunks.append(line.strip()) else: # v4: host mode filters a Caudoviricetes download (older "hosts_"/"hosts3_" sets missed many phages) tag = "all_Caudoviricetes" if mode == "all" else "hosts4_" + "_".join(sorted(h.lower() for h in hosts))[:120] stem = os.path.join(ncbi_dir, f"ncbi_phages_{re.sub(r'[^A-Za-z0-9_.-]+', '_', tag)}") age = ncbi_age_days(stem) if os.path.isfile(stem + ".meta.tsv") else None if age is not None and age > NCBI_MAX_AGE_DAYS and not refresh: print(f"The NCBI genomes were downloaded {age:.0f} days ago; downloading them again to include new phages…", flush=True) refresh = True if refresh or not os.path.isfile(stem + ".meta.tsv"): print("Downloading complete phage genomes from NCBI (once; refresh with the generator's checkbox)…", flush=True) ncbi_download(ncbi_dir, mode, hosts, stem) for ext in (".asm20.mmi",): if os.path.exists(stem + ext): os.remove(stem + ext) index = stem + ".asm20.mmi" if not os.path.isfile(index): print("Indexing the NCBI phage genomes for minimap2 (once per download)…", flush=True) r = subprocess.run(["minimap2", "-x", "asm20", "-t", str(threads), "-d", index + ".part", stem + ".fasta"], capture_output=True, text=True) if r.returncode != 0: raise RuntimeError("minimap2 could not index the NCBI genomes:\\n" + r.stderr[-800:]) os.replace(index + ".part", index) meta = {} with open(stem + ".meta.tsv") as fh: next(fh) for line in fh: f = line.rstrip("\\n").split("\\t") if f and f[0]: desc = f[1] if len(f) > 1 else "" host = f[2] if len(f) > 2 else "" if not host: # e.g. "Staphylococcus phage X, complete genome" -> Staphylococcus mh = re.match(r"\\s*([A-Z][a-z]+)\\s+(?:sp\\.\\s+)?phage\\b", desc) host = mh.group(1) if mh else "" meta[f[0]] = {"description": desc, "host": host, "length": f[3] if len(f) > 3 else ""} info = {} if os.path.isfile(stem + ".info.txt"): for line in open(stem + ".info.txt"): k, _, v = line.rstrip("\\n").partition("\\t") info[k] = v return index, meta, info def paf_lines(cmd, what): """Run minimap2 and return only the 12 standard PAF columns, read line by line. v2.1: the CIGAR strings of -c alignments of hundreds of genomes can take several GB when the whole output is held in memory, which can stop WSL.""" import tempfile out = [] with tempfile.TemporaryFile(mode="w+") as err: p = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=err, text=True, bufsize=1 << 20) for line in p.stdout: f = line.split("\\t", 12)[:12] if len(f) == 12: out.append(f) if p.wait() != 0: err.seek(0) raise RuntimeError(f"minimap2 ({what}) failed:\\n" + err.read()[-800:]) return out def mem_note(): """free memory, for the progress lines (helps when WSL runs out of memory)""" try: for line in open("/proc/meminfo"): if line.startswith("MemAvailable:"): return f" [{int(line.split()[1]) / 1048576:.1f} GB free]" except OSError: pass return "" def gpos(qoff, name): """contig query name -> (genome id, offset of the contig in the genome's coordinates)""" return qoff.get(name, (name, 0)) if qoff else (name, 0) def search(index, qpath, threads, qoff=None): """minimap2 (asm20) of all query contigs against an index -> {genome: {target: totals}} (contigs of one genome are placed one after the other, see genome_query).""" hits = {} for f in paf_lines(["minimap2", "-c", "-x", "asm20", "-N", "100", "-p", "0.2", "-t", str(threads), index, qpath], "search"): g, off = gpos(qoff, f[0]) d = hits.setdefault(g, {}).setdefault(f[5], {"tlen": int(f[6]), "iv": [], "tiv": [], "st": [], "match": 0, "block": 0}) d["st"].append(f[4]) d["iv"].append((int(f[2]) + off, int(f[3]) + off)) d["tiv"].append((int(f[7]), int(f[8]))) d["match"] += int(f[9]) d["block"] += int(f[10]) return hits def rank_hits(hits_q, L): ranked = [] for tid, d in (hits_q or {}).items(): cov = sum(e - s for s, e in merge(d["iv"])) tcov = sum(e - s for s, e in merge(d["tiv"])) ident = 100.0 * d["match"] / d["block"] if d["block"] else 0 ranked.append({"score": cov * ident, "id": tid, "tlen": d["tlen"], "ident": round(ident, 2), "qcov": round(100.0 * cov / L, 1) if L else 0, "aln": cov, "taln": tcov, "scov": round(100.0 * tcov / d["tlen"], 1) if d["tlen"] else 0}) ranked.sort(key=lambda x: -x["score"]) return ranked def contig_hits(hq, offs, meta): """v3.0: the closest PhagesDB phage of each searched contig: [[contig, length, phage, cluster, % of the contig aligned, identity]].""" out = [] for a, b, name, ln in offs: best = None for tid, d in (hq or {}).items(): iv = [(max(s, a), min(e, b)) for s, e in d["iv"] if e > a and s < b] if not iv: continue cov = sum(e - s for s, e in merge(iv)) ident = 100.0 * d["match"] / d["block"] if d["block"] else 0 if not best or cov * ident > best[0]: best = (cov * ident, tid, cov, ident) if best and best[2] >= min(1000, 0.2 * ln): m = meta.get(best[1], {}) out.append([name, ln, m.get("phage") or best[1], m.get("cluster", ""), round(100.0 * best[2] / ln, 1), round(best[3], 2)]) else: out.append([name, ln, "", "", 0, 0]) return out def contig_of(offs, pos): return next((i for i, o in enumerate(offs) if o[0] <= pos < o[1]), 0) def main_contigs(d, offs, tlen): """v3.1: the contigs that make up the phage genome matching the closest relative (the rest of the sample's contigs are host DNA, debris or another phage). Returns (indexes, % of the relative covered, % of those contigs aligned, their length, aligned bases, relative bases covered).""" per = {} for (qs, qe), (ts, te) in zip(d["iv"], d["tiv"]): k = contig_of(offs, qs) p = per.setdefault(k, {"q": [], "t": []}) p["q"].append((qs, qe)) p["t"].append((min(ts, te), max(ts, te))) cov = lambda iv: sum(e - s for s, e in merge(iv)) order = sorted(per, key=lambda k: -cov(per[k]["t"])) chosen, tiv = [], [] for k in order: ln = offs[k][1] - offs[k][0] gain = cov(tiv + per[k]["t"]) - cov(tiv) if not chosen or (gain >= max(1000, 0.02 * tlen) and cov(per[k]["q"]) >= 0.3 * ln): chosen.append(k) tiv += per[k]["t"] if not chosen: return [], 0.0, 0.0, 0, 0, 0 Lm = sum(offs[k][1] - offs[k][0] for k in chosen) qa = sum(cov(per[k]["q"]) for k in chosen) ta = cov(tiv) return sorted(chosen), round(100.0 * ta / tlen, 1) if tlen else 0.0, round(100.0 * qa / Lm, 1) if Lm else 0.0, Lm, qa, ta def cluster_letter(c): c = str(c or "") if c.lower() == "singleton": return "Singleton" m = re.match(r"[A-Z]+", c) return m.group(0) if m else "" def cluster_call(hq, offs, chosen, Lm, meta, clu_len): """v3.1: the PhagesDB cluster from the hits on this genome's own contigs only (not on host DNA or another sample's phage), with how many of the top hits agree and the SEA-PHAGES rule (>= 50% span to a cluster member).""" rng = [(offs[k][0], offs[k][1]) for k in chosen] ranked = [] for tid, d in (hq or {}).items(): q, t = [], [] for (qs, qe), (ts, te) in zip(d["iv"], d["tiv"]): if any(a <= qs < b for a, b in rng): q.append((qs, qe)) t.append((min(ts, te), max(ts, te))) if not q: continue qa = sum(e - s for s, e in merge(q)) ta = sum(e - s for s, e in merge(t)) ident = 100.0 * d["match"] / d["block"] if d["block"] else 0 ranked.append((qa * ident, tid, qa, ta, d["tlen"], ident)) ranked.sort(key=lambda x: -x[0]) out = {"cluster": "", "letter": "", "agree": "", "span": 0.0, "note": "", "hit": "", "ident": "", "range": ""} if not ranked or not Lm or ranked[0][2] < 0.05 * Lm: out["note"] = "no PhagesDB phage matches this genome's contigs: not an Actinobacteriophage, or a new kind" return out top = ranked[:5] m0 = meta.get(top[0][1], {}) c0 = m0.get("cluster", "") or "?" let = cluster_letter(c0) span = max(100.0 * top[0][2] / Lm, 100.0 * top[0][3] / top[0][4] if top[0][4] else 0) same = [t for t in top if cluster_letter(meta.get(t[1], {}).get("cluster", "")) == let] subs = sorted({meta.get(t[1], {}).get("cluster", "") for t in same} - {""}) out.update({"letter": let, "agree": f"{len(same)} of {len(top)}", "span": round(span, 1), "hit": m0.get("phage") or top[0][1], "ident": round(top[0][5], 2)}) lens = [x for s in subs for x in clu_len.get(s, [])] if len(lens) >= 3: out["range"] = [min(lens), max(lens)] disp = c0 if len(subs) <= 1 else f"{let} ({' or '.join(subs)})" if span >= 50: out["cluster"] = disp out["note"] = (f"{len(same)} of the {len(top)} closest PhagesDB phages are cluster {let}; {round(span)}% span to " f"{out['hit']} ({c0}), over the 50% needed to place it") else: out["cluster"] = disp + "?" out["note"] = (f"only {round(span)}% span to {out['hit']} ({c0}), under the 50% needed to place it in a cluster: " f"possibly a new subcluster or cluster (check gene content similarity, 35% or more)") return out def gb_group(desc): """'Staphylococcus phages (GenBank)' from 'Staphylococcus phage 80B' or 'Staphylococcus virus BESEP5'.""" m = re.match(r"\\s*([A-Z][a-z]+)\\s+(phage|virus)\\b", str(desc or "")) return f"{m.group(1)} phages (GenBank)" if m else "" def completeness(n_main, circ, Lm, rng, scov, cvq, let): """v3.1: completeness from several lines of evidence (CheckV alone underrates new Actinobacteriophages).""" ev = [["one contig", n_main == 1], ["circular contig", True if circ == "yes" else (False if circ == "no" else None)], [f"length within cluster {let}'s range ({rng[0] / 1000:.1f}–{rng[1] / 1000:.1f} kb)" if rng else "length in its cluster's range", (0.97 * rng[0] <= Lm <= 1.03 * rng[1]) if rng else None], ["covers 90% or more of the closest phage", (float(scov) >= 90) if scov != "" else None], ["CheckV high-quality", (cvq in ("Complete", "High-quality")) if cvq else None]] ok = sum(1 for e in ev if e[1]) known = sum(1 for e in ev if e[1] is not None) if known < 2: v = "not enough evidence" elif n_main == 1 and (ev[1][1] or ev[2][1] or ev[3][1]): v = "complete" elif ok >= 3: v = "likely complete" elif n_main > 1 and ev[3][1]: v = "complete, in pieces" else: v = "incomplete" return v, ev, f"{ok} of {known}" def _group(r): if r.get("cluster_letter"): return "PDB:" + r["cluster_letter"] return ("GB:" + r["gb_group"]) if r.get("gb_group") else "" def post_rows(rows, checks): """v3.1: (1) 'same phage' pairs whose genomes belong to different clusters are DNA of one sample's phage inside the other sample (barcode bleed or cross-contamination), not the same phage; (2) Ligase / No-Ligase pairs: which assembly is better; then the notes and the run log.""" by = {r["sample"]: r for r in rows} for r in rows: keep = [] for n, v in r["dups"]: o = by.get(n) ga, gb = _group(r), _group(o) if o else "" if o and ga and gb and ga != gb: if (r["other_contigs"][1], r["searched"]) > (o["other_contigs"][1], o["searched"]): lab = o["cluster_letter"] and "cluster " + o["cluster_letter"] or o.get("gb_group", "") r["contains"].append([n, v, lab]) continue keep.append([n, v]) r["dups"] = keep rank = {"complete": 0, "likely complete": 1, "complete, in pieces": 2, "incomplete": 3} pairs = {} for r in rows: k = re.sub(r"[-_ ]?(no[-_ ]?ligase|ligase)$", "", r["sample"], flags=re.I) if k != r["sample"]: pairs.setdefault(k.lower(), []).append(r) for k, lst in pairs.items(): if len(lst) != 2: continue def key(r): b = r["best"] or {} return (rank.get(r["completeness"], 9), len(r["main_contigs"]), r["other_contigs"][1], abs(r["main_len"] - (b.get("slen") or r["main_len"])), -float(b.get("scov") or 0)) a, b2 = sorted(lst, key=key) if key(a) == key(b2): continue why = [] if len(a["main_contigs"]) != len(b2["main_contigs"]): why.append(f"{len(a['main_contigs'])} contig vs {len(b2['main_contigs'])}") if a["completeness"] != b2["completeness"]: why.append(f"{a['completeness']} vs {b2['completeness']}") if a["other_contigs"][1] != b2["other_contigs"][1] and not why: why.append(f"{a['other_contigs'][1]:,} vs {b2['other_contigs'][1]:,} bp of other contigs") if not why: why.append("closer in length to the closest phage") a["pair"] = {"with": b2["sample"], "better": True, "why": "; ".join(why)} b2["pair"] = {"with": a["sample"], "better": False, "why": "; ".join(why)} for r in rows: flags = [] if r["host_pct"] != "" and float(r["host_pct"]) >= 50: flags.append(f"mostly host DNA ({r['host_pct']}% matches {r['host_name']})") elif r["host_pct"] != "" and float(r["host_pct"]) >= 5: flags.append(f"{r['host_pct']}% host DNA") if r["cv_quality"] in ("Low-quality", "Not-determined"): flags.append(f"CheckV: {r['cv_quality']}") if r["dups"]: flags.append("same phage as " + ", ".join(f"{n} ({v}%)" for n, v in r["dups"][:3])) if r["contains"]: flags.append("contains DNA of " + ", ".join(f"{n} ({lab})" for n, v, lab in r["contains"][:3]) + ": barcode bleed or cross-contamination") if r["pair"] and not r["pair"]["better"]: flags.append(f"{r['pair']['with']} is the better assembly ({r['pair']['why']})") if r["host_mismatch"]: flags.append(r["host_mismatch"]) if flags: r["note"] = "; ".join([r["note"]] + flags if r["note"] else flags) bst = r["best"] status = (f"{bst['name']} [{bst['db']}] {bst['ident']}% identity, {bst['qcov']}% of yours covered, " f"covers {bst['scov']}% of it ({bst['slen']:,} bp) — {r['status']}, {bst['rank']} (~{bst['sim']}%)" if bst else "none") print(f" {r['sample']:<28} {status}{' — ' + r['note'] if r['note'] else ''}") if r["searched"] > 1: print(f" {'':<28} {r['searched']} contigs searched; this genome: {len(r['main_contigs'])} contig(s), {r['main_len']:,} bp" + (f"; {r['other_contigs'][0]} other contig(s), {r['other_contigs'][1]:,} bp" if r["other_contigs"][0] else "")) cc = r["cluster_call"] print(f" {'':<28} cluster: {cc['cluster'] or r['gb_group'] or '—'} ({cc['agree'] + ' top hits agree' if cc['agree'] else cc['note']}); " f"completeness: {r['completeness']} ({r['completeness_score']} lines of evidence)") if r["_ties_n"]: print(f" {'':<28} {r['_ties_n']} other phage(s) within 0.5% of the best: " + ", ".join(t[0] for t in r["ties"][:4])) if checks["host"] or checks["checkv"]: bits = [] if checks["host"]: bits.append(f"host DNA {r['host_pct']}%") if r["cv_quality"]: bits.append(f"CheckV {r['cv_quality']}, completeness {r['cv_completeness']}%") print(f" {'':<28} " + "; ".join(bits)) if r["ncbi"]: nb = r["ncbi"] print(f" {'':<28} GenBank: {nb['acc']} {nb['desc'][:60]} — {nb['ident']}% ident, {nb['qcov']}% query cover") elif r["ncbi"] == {}: print(f" {'':<28} GenBank: no phage aligns") r.pop("_ties_n", None) def aln_blocks(d, offs, min_len=150, cap=400): """Where each searched contig aligns on the relative: [[start, end, contig index, strand], ...] (merged).""" by = {} for (qs, qe), (ts, te), st in zip(d["iv"], d["tiv"], d.get("st") or ["+"] * len(d["iv"])): k = next((i for i, o in enumerate(offs) if o[0] <= qs < o[1]), 0) by.setdefault((k, st), []).append((ts, te)) out = [[a, b, k, st] for (k, st), iv in by.items() for a, b in merge(iv) if b - a >= min_len] out.sort(key=lambda x: (x[0], x[1])) return out[:cap] def ictv_sim(ident, aln_q, aln_t, la, lb): """Approximate VIRIDIC-style similarity: 2 x identical bases / (both lengths) x 100.""" if not la + lb: return 0.0 return round(min(100.0, 200.0 * (float(ident) / 100.0) * min(aln_q, aln_t) / (la + lb)), 1) def ictv_rank(v): return "same species" if v >= 95 else ("same genus" if v >= 70 else "different genus") def url_ncbi(acc): return "https://www.ncbi.nlm.nih.gov/nuccore/" + urllib.parse.quote(str(acc)) # ---------------------------------------------------------------- genome checks def host_check(qpath, host_paths, threads, qoff=None, glen=None): """Share of each query genome aligned (>= 85% identity) to the host genome(s) -> {qid: (pct, host label)}.""" import tempfile labels, out = [], {} with tempfile.TemporaryDirectory() as tmp: hf = os.path.join(tmp, "hosts.fa") with open(hf, "w") as fo: for i, hp in enumerate(host_paths, 1): label = re.sub(r"\\.(fasta|fa|fna|fas)(\\.gz)?$", "", os.path.basename(hp), flags=re.I) labels.append(label) for name, _, seq in read_fasta(hp): write_fasta(fo, f"H{i}__{name}", seq) paf = paf_lines(["minimap2", "-c", "-x", "asm20", "-t", str(threads), hf, qpath], "host check") spans, per_host, qlen = {}, {}, {} for f in paf: if int(f[10]) == 0 or int(f[9]) / int(f[10]) < 0.85: continue h = int(f[5].split("__", 1)[0][1:]) - 1 q, off = gpos(qoff, f[0]) qlen[q] = (glen or {}).get(q) or int(f[1]) spans.setdefault(q, []).append((int(f[2]) + off, int(f[3]) + off)) per_host.setdefault(q, {}).setdefault(h, []).append((int(f[2]) + off, int(f[3]) + off)) for q, iv in spans.items(): cov = sum(e - s for s, e in merge(iv)) best = max(per_host[q].items(), key=lambda kv: sum(e - s for s, e in merge(kv[1])))[0] out[q] = (round(100.0 * cov / qlen[q], 1), labels[best]) return out def find_checkv_db(db): """The CheckV database folder (the one holding genome_db and hmm_db), at DB or one level below.""" if not db or not os.path.isdir(db): return "" if os.path.isdir(os.path.join(db, "genome_db")) and os.path.isdir(os.path.join(db, "hmm_db")): return db for sub in sorted(os.listdir(db), reverse=True): if sub.endswith(".incomplete"): continue d = os.path.join(db, sub) if os.path.isdir(os.path.join(d, "genome_db")) and os.path.isdir(os.path.join(d, "hmm_db")): return d return "" def run_checkv(qpath, out_dir, bin_dir, db, download, threads): """CheckV end_to_end on the query genomes -> {qid: quality_summary row}.""" env = dict(os.environ) if bin_dir: env["PATH"] = bin_dir + os.pathsep + env.get("PATH", "") import shutil if not shutil.which("checkv", path=env["PATH"]): raise RuntimeError("CheckV was not found" + (f" in {bin_dir}" if bin_dir else "")) dbdir = find_checkv_db(db) if not dbdir: if not download: raise RuntimeError(f"no CheckV database in {db}") os.makedirs(db, exist_ok=True) print(f"Downloading the CheckV database into {db} (once; a few GB)…", flush=True) r = subprocess.run(["checkv", "download_database", db], env=env, capture_output=True, text=True) dbdir = find_checkv_db(db) if r.returncode != 0 or not dbdir: raise RuntimeError("CheckV database download failed:\\n" + (r.stderr or r.stdout)[-800:]) # v2.3: the DIAMOND file is built at the end of "checkv download_database"; an interrupted download (or a manual # one) leaves it missing -> build it from checkv_reps.faa, or download the database again gdb = os.path.join(dbdir, "genome_db") dmnd = os.path.join(gdb, "checkv_reps.dmnd") if not os.path.isfile(dmnd) or os.path.getsize(dmnd) == 0: faa = next((os.path.join(gdb, f) for f in ("checkv_reps.faa", "checkv_reps.faa.gz") if os.path.isfile(os.path.join(gdb, f))), "") if faa and shutil.which("diamond", path=env["PATH"]): print(f"The CheckV database is missing {os.path.basename(dmnd)}; building it with DIAMOND (once; a few minutes)…", flush=True) r = subprocess.run(["diamond", "makedb", "--in", faa, "--db", os.path.join(gdb, "checkv_reps"), "--threads", str(threads)], env=env, capture_output=True, text=True) if r.returncode != 0 or not os.path.isfile(dmnd): raise RuntimeError("could not build the CheckV DIAMOND database:\\n" + (r.stderr or r.stdout)[-600:] + f"\\nBuild it by hand: cd {gdb} && diamond makedb --in {os.path.basename(faa)} --db checkv_reps") elif download: print(f"The CheckV database in {dbdir} is incomplete; downloading it again (once; a few GB)…", flush=True) bad = dbdir.rstrip("/") + ".incomplete" if os.path.exists(bad): shutil.rmtree(bad, ignore_errors=True) os.replace(dbdir, bad) r = subprocess.run(["checkv", "download_database", db], env=env, capture_output=True, text=True) dbdir = find_checkv_db(db) if r.returncode != 0 or not dbdir or not os.path.isfile(os.path.join(dbdir, "genome_db", "checkv_reps.dmnd")): raise RuntimeError("CheckV database download failed:\\n" + (r.stderr or r.stdout)[-800:]) shutil.rmtree(bad, ignore_errors=True) else: raise RuntimeError(f"the CheckV database in {dbdir} is incomplete ({os.path.basename(dmnd)} and checkv_reps.faa are " "missing); tick 'Download the database if it is missing' or delete the folder and download it again") work = os.path.join(out_dir, "checkv") print(f"Running CheckV (database {os.path.basename(dbdir)})…", flush=True) r = subprocess.run(["checkv", "end_to_end", qpath, work, "-d", dbdir, "-t", str(threads)], env=env, capture_output=True, text=True) summ = os.path.join(work, "quality_summary.tsv") if r.returncode != 0 or not os.path.isfile(summ): raise RuntimeError("CheckV failed:\\n" + (r.stderr or r.stdout)[-800:]) res = {} with open(summ) as fh: cols = next(fh).rstrip("\\n").split("\\t") for line in fh: d = dict(zip(cols, line.rstrip("\\n").split("\\t"))) res[d.get("contig_id", "")] = d return res def similarity_matrix(qpath, info, threads, qoff=None): """VIRIDIC-style intergenomic similarity between all pairs of genomes: SIM = 2 x identical bases in the alignments / (length A + length B) x 100, from minimap2 (asm20) alignments (overlapping alignments counted once).""" alns = {} for f in paf_lines(["minimap2", "-c", "-x", "asm20", "-X", "-t", str(threads), qpath, qpath], "similarity"): (g1, o1), (g2, o2) = gpos(qoff, f[0]), gpos(qoff, f[5]) if g1 == g2: continue a, b = sorted((g1, g2)) qs, qe, ts, te = int(f[2]) + o1, int(f[3]) + o1, int(f[7]) + o2, int(f[8]) + o2 if a != g1: qs, qe, ts, te = ts, te, qs, qe alns.setdefault((a, b), []).append((int(f[9]), qs, qe, ts, te, int(f[10]))) sim, pair = {}, {} for (a, b), lst in alns.items(): lst.sort(key=lambda x: -x[0]) used_a, used_b, total, blk = [], [], 0, 0 for m, qs, qe, ts, te, bl in lst: ov_a = sum(max(0, min(qe, e) - max(qs, s)) for s, e in used_a) ov_b = sum(max(0, min(te, e) - max(ts, s)) for s, e in used_b) if ov_a > 0.5 * (qe - qs) or ov_b > 0.5 * (te - ts): continue used_a.append((qs, qe)) used_b.append((ts, te)) total += m blk += bl la, lb = info[a]["length"], info[b]["length"] sim[(a, b)] = round(min(100.0, 200.0 * total / (la + lb)), 1) if la + lb else 0.0 ca, cb = sum(e - s for s, e in merge(used_a)), sum(e - s for s, e in merge(used_b)) short = 100.0 * (ca / la if la <= lb else cb / lb) if la and lb else 0.0 pair[(a, b)] = (round(100.0 * total / blk, 2) if blk else 0.0, round(min(100.0, short), 1)) return sim, pair def groups_at(qids, sim, cut): """Connected groups of genomes with similarity >= cut (single linkage).""" parent = {q: q for q in qids} def find(x): while parent[x] != x: parent[x] = parent[parent[x]] x = parent[x] return x for (a, b), s in sim.items(): if s >= cut: parent[find(a)] = find(b) comp = {} for q in qids: comp.setdefault(find(q), []).append(q) return [sorted(v) for v in comp.values()] # ---------------------------------------------------------------- input genomes # folders in a pipeline output folder that are not phages (reports, uploads, indexes, earlier runs of this tool) NOT_PHAGE_DIR = re.compile(r"^(\\.|closest_phages|assembly_comparison|nanocomp_|sharepoint_upload|phagesdb_index|" r"sample_logs|timings_|upload)", re.I) def pick_fasta(folder, name, pattern): files = sorted(f for f in os.listdir(folder) if f.lower().endswith(FASTA_EXT) and not f.startswith(".") and os.path.isfile(os.path.join(folder, f))) if not files: return None if pattern: rx = re.compile("^" + re.escape(pattern).replace(r"\\*", ".*").replace(r"\\{phage\\}", re.escape(name)) + "$", re.I) hit = [f for f in files if rx.match(f)] return os.path.join(folder, hit[0]) if hit else None low = {f.lower(): f for f in files} for cand in (f"{name}.fasta", f"{name}.fa", f"{name}.fna", f"{name}.longest_contig.fasta"): if cand.lower() in low: return os.path.join(folder, low[cand.lower()]) skip = re.compile(r"(\\.flye|\\.autocycler|unrotated|_reoriented|\\.ids)\\.f", re.I) rest = [f for f in files if not skip.search(f)] or files return os.path.join(folder, max(rest, key=lambda f: os.path.getsize(os.path.join(folder, f)))) NO_FASTA = [] def find_genomes(root, pattern, skip_dirs=()): out, skipped = [], [] del NO_FASTA[:] skip = {os.path.realpath(d) for d in skip_dirs if d} for entry in sorted(os.listdir(root), key=str.lower): p = os.path.join(root, entry) if os.path.isdir(p) and (NOT_PHAGE_DIR.match(entry) or os.path.realpath(p) in skip): if not entry.startswith("."): skipped.append(entry) continue if os.path.isdir(p): f = pick_fasta(p, entry, pattern) if f: out.append((entry, f)) else: NO_FASTA.append(entry) elif os.path.isfile(p) and entry.lower().endswith(FASTA_EXT) and not pattern and not entry.startswith("."): stem = re.sub(r"\\.(fasta|fa|fna|fas)(\\.gz)?$", "", entry, flags=re.I) out.append((stem, p)) return out, skipped def genome_query(name, path, recs, min_len): """One query per genome: its contigs >= min_len joined with N spacers (the longest contig if none is that long).""" recs = sorted(recs, key=lambda r: -len(r[2])) big = [r for r in recs if len(r[2]) > MAX_CONTIG] use = [r for r in recs if min_len <= len(r[2]) <= MAX_CONTIG] or [r for r in recs if len(r[2]) <= MAX_CONTIG][:1] or recs[:1] offs, parts, pos = [], [], 0 for r in use: if parts: parts.append(SPACER) pos += len(SPACER) offs.append((pos, pos + len(r[2]), r[0], len(r[2]))) parts.append(r[2]) pos += len(r[2]) L = sum(o[3] for o in offs) gcn = sum(r[2].upper().count("G") + r[2].upper().count("C") for r in use) if big: print(f" {name}: {len(big)} contig(s) over {MAX_CONTIG // 1000:,} kb not searched (host DNA, not phage)") return {"sample": name, "file": path, "contigs": len(recs), "searched": len(use), "length": L, "longest": len(recs[0][2]), "circ_hdr": header_topology(recs[0][1]), "offs": offs, "gc": round(100.0 * gcn / L, 1) if L else 0, "seq": "".join(parts)} def header_topology(desc): """'yes'/'no' from a FASTA header tag (circular=true, Autocycler/Flye style), else ''.""" m = re.search(r"\\bcirc(?:ular)?\\s*[=:]\\s*(true|yes|y|false|no|n)\\b", str(desc), re.I) return "" if not m else ("yes" if m.group(1).lower() in ("true", "yes", "y") else "no") def norm_key(name): return re.sub(r"[^a-z0-9]+", "", str(name).lower()) def read_isolation_hosts(paths): """phage -> isolation host from pasted lists or CSV/TSV exports (e.g. PhageTracker): 'PheterStaphylococcus epidermidis', 'Pheter, S. epidermidis' or a table with phage/name and host columns.""" import csv out = {} for p in paths: if not p: continue if not os.path.isfile(p): print(f"⚠️ Isolation host list not found: {p}") continue lines = [ln for ln in open_text(p).read().splitlines() if ln.strip() and not ln.lstrip().startswith("#")] if not lines: continue def split(ln): # separator per line: tab, comma, semicolon, colon, else two spaces / "name Genus species" d = next((c for c in ("\\t", ",", ";", ":") if c in ln), None) if d: return next(csv.reader([ln], delimiter=d)) return re.split(r"\\s{2,}|\\s+(?=[A-Z][a-z]+\\s)", ln.strip(), maxsplit=1) table = [split(ln) for ln in lines] head = [c.strip().lower() for c in table[0]] pc = hc = None for i, c in enumerate(head): if pc is None and re.fullmatch(r"(phage|phage[ _]?name|name|sample|phage[ _]?id)", c): pc = i if hc is None and re.search(r"isolation[ _]?host|^host([ _]?(species|strain|name|organism|genus))?$", c): hc = i if pc is not None and hc is not None: body = table[1:] else: pc, hc, body = 0, 1, table for r in body: if len(r) > max(pc, hc) and r[pc].strip() and r[hc].strip(): out[norm_key(r[pc])] = r[hc].strip() return out def host_genus(h): h = str(h or "").strip() if not h: return "" w = re.sub(r"[^A-Za-z-]", "", h.split()[0]) w = GENUS_SYNONYMS.get(w.lower(), w) return w if len(w) > 2 else h # genera split from Mycobacterium in 2018 (PhagesDB lists M. smegmatis as Mycolicibacterium smegmatis) GENUS_SYNONYMS = {"mycolicibacterium": "Mycobacterium", "mycobacteroides": "Mycobacterium", "mycolicibacter": "Mycobacterium", "mycolicibacillus": "Mycobacterium"} def read_kv(path): d = {} if os.path.isfile(path): for line in open(path, errors="replace"): k, _, v = line.rstrip("\\n").partition("\\t") d[k] = v return d def pipeline_info(name, fasta, out_dir, in_dir=""): """coverage.html link (relative to the report, and relative to the input folder for a SharePoint base URL), warnings and circularity (Flye/Autocycler topology of the longest contig) from a Nanopore phage pipeline output folder.""" d = os.path.dirname(os.path.abspath(fasta)) link, rel, flags = "", "", [] for cand in (f"{name}.coverage.html", os.path.basename(d) + ".coverage.html"): p = os.path.join(d, cand) if os.path.isfile(p): link = os.path.relpath(p, os.path.abspath(out_dir)).replace(os.sep, "/") rel = os.path.relpath(p, os.path.abspath(in_dir)).replace(os.sep, "/") if in_dir else cand break circ = {"Y": "yes", "N": "no"}.get(read_kv(os.path.join(d, "qc", "longest_contig.tsv")).get("flye_circular", "").strip(), "") f1 = read_kv(os.path.join(d, "qc", "sample_flags.txt")).get("flags", "") for x in re.split(r";\\s*", f1): if x.strip() and x.strip().lower() not in ("none", "ok", "-"): flags.append(x.strip()) f2 = read_kv(os.path.join(d, "qc", "input_check.txt")).get("flag", "") if f2.strip() and f2.strip().lower() != "none": flags.append("input check: " + f2.strip()) return link, flags, circ, rel def pipeline_reads(fasta): """Bases sequenced (deduplicated; host reads removed when that was done) and read depth along the genome, from a Nanopore phage pipeline output folder: qc/readstats_*.json, qc/coverage_stats.txt, qc/coverage_dips.tsv.""" qc = os.path.join(os.path.dirname(os.path.abspath(fasta)), "qc") out = {"bases": "", "reads_label": "", "med_depth": "", "pct_low": "", "dips": ""} for fn, label in (("readstats_host_depleted.json", "phage reads (host reads removed)"), ("readstats_raw.json", "deduplicated reads")): p = os.path.join(qc, fn) if os.path.isfile(p): try: tb = json.load(open(p)).get("total_bases") if tb: out["bases"], out["reads_label"] = int(tb), label break except (ValueError, OSError): pass cs = read_kv(os.path.join(qc, "coverage_stats.txt")) out["med_depth"], out["pct_low"] = cs.get("median_depth", ""), cs.get("pct_below_lowcov", "") out["dips"] = read_kv(os.path.join(qc, "coverage_dips.tsv")).get("regions", "") return out def sequenced(folder): """Reads and bases sequenced for a phage (all its input FASTQs), from a pipeline phage folder: qc/read_sets.tsv (written for assembled and failed phages), else qc/readstats_raw.json (deduplicated reads). Also the failed step and the failure report of a phage that was not assembled.""" qc = os.path.join(folder, "qc") out = {"reads": None, "bases": None, "n50": "", "mean_q": "", "src": "", "step": "", "report": ""} rs = os.path.join(qc, "read_sets.tsv") if os.path.isfile(rs): try: with open(rs, errors="replace") as fh: head = fh.readline().rstrip("\\n").split("\\t") rows = [dict(zip(head, ln.rstrip("\\n").split("\\t"))) for ln in fh if ln.strip()] allr = [r for r in rows if r.get("role") == "all_inputs"] ins = [r for r in rows if r.get("role") == "input"] use = allr[:1] or ins if use: out["reads"] = sum(int(r.get("reads") or 0) for r in use) out["bases"] = sum(int(r.get("bases") or 0) for r in use) if len(use) == 1: out["n50"], out["mean_q"] = use[0].get("n50", ""), use[0].get("mean_q", "") else: qw = [(float(r["mean_q"]), int(r.get("reads") or 0)) for r in use if r.get("mean_q")] if qw and sum(w for _, w in qw): out["mean_q"] = round(sum(q * w for q, w in qw) / sum(w for _, w in qw), 2) out["src"] = "all input reads" except (OSError, ValueError): pass if out["reads"] is None: p = os.path.join(qc, "readstats_raw.json") if os.path.isfile(p): try: j = json.load(open(p)) out["reads"], out["bases"] = j.get("read_count"), j.get("total_bases") out["n50"], out["mean_q"] = j.get("read_n50") or "", j.get("mean_q_score") or "" out["src"] = "deduplicated reads" except (OSError, ValueError): pass out["step"] = read_kv(os.path.join(qc, "fail_report.cfg")).get("failed_step", "") name = os.path.basename(os.path.normpath(folder)) if os.path.isfile(os.path.join(folder, f"{name}.coverage.html")): out["report"] = f"{name}/{name}.coverage.html" return out def _ranks(v): o = sorted(range(len(v)), key=lambda i: v[i]) r = [0.0] * len(v) i = 0 while i < len(o): j = i while j + 1 < len(o) and v[o[j + 1]] == v[o[i]]: j += 1 for k in range(i, j + 1): r[o[k]] = (i + j) / 2.0 i = j + 1 return r def spearman(x, y): rx, ry = _ranks(x), _ranks(y) n = len(x) m = (n - 1) / 2.0 sxy = sum((a - m) * (b - m) for a, b in zip(rx, ry)) sxx = sum((a - m) ** 2 for a in rx) syy = sum((b - m) ** 2 for b in ry) return sxy / (sxx * syy) ** 0.5 if sxx and syy else 0.0 def gc_bias(fasta, sample, win=500, min_contig=5000): """v2.4: read depth vs GC content in 500 bp windows, from the pipeline's per-base coverage table (qc/.coverage.tsv.gz: contig, pos, consensus base, depth). Contigs of 5 kb or more; the first and last window of each contig are left out (ends have lower depth anyway). Depth is relative to the phage's median. Returns rho (Spearman, GC vs depth), slope (log2 depth change per +10% GC), and per-phage medians in 1% bins of GC minus the phage's mean GC, so phages with different GC content can be compared.""" import math qc = os.path.join(os.path.dirname(os.path.abspath(fasta)), "qc") p = os.path.join(qc, f"{sample}.coverage.tsv.gz") if not os.path.isfile(p): cand = [f for f in os.listdir(qc) if f.endswith(".coverage.tsv.gz")] if os.path.isdir(qc) else [] if len(cand) != 1: return None p = os.path.join(qc, cand[0]) wins, cur, acc = [], None, [] def flush(): if len(acc) >= max(3, min_contig // win): wins.extend(acc[1:-1]) gc = dep = n = 0 try: with gzip.open(p, "rt") as fh: next(fh, None) for line in fh: f = line.split("\\t", 5) if len(f) < 4: continue if f[0] != cur: flush() cur, acc, gc, dep, n = f[0], [], 0, 0, 0 b = f[2][:1].upper() gc += b in "GC" dep += int(f[3]) n += 1 if n == win: acc.append((100.0 * gc / win, dep / win)) gc = dep = n = 0 flush() except (OSError, ValueError, EOFError): return None if len(wins) < 20: return None ds = sorted(w[1] for w in wins) med = ds[len(ds) // 2] if med < 10: return None x = [w[0] for w in wins] y = [w[1] / med for w in wins] rho = spearman(x, y) ly = [math.log2(max(v, 0.01)) for v in y] mx, my = sum(x) / len(x), sum(ly) / len(ly) sxx = sum((a - mx) ** 2 for a in x) slope = 10 * sum((a - mx) * (b - my) for a, b in zip(x, ly)) / sxx if sxx else 0.0 bins = {} for a, v in zip(x, y): k = max(-15, min(15, int(round(a - mx)))) bins.setdefault(k, []).append(v) bout = {str(k): [round(sorted(v)[len(v) // 2], 3), len(v)] for k, v in bins.items() if len(v) >= 3} low = [v for a, v in zip(x, y) if a <= mx - 5] return {"rho": round(rho, 3), "slope": round(slope, 3), "n": len(wins), "gc": round(mx, 1), "med": round(med, 1), "bins": bout, "low5": round(sorted(low)[len(low) // 2], 3) if len(low) >= 5 else None} def check_mode(args, genomes, skipped): """--check: list the genomes and FASTA files that would be searched, without searching.""" iso = read_isolation_hosts(args.isolation_hosts) print(f"\\n{len(genomes)} genome(s) in {args.input} (contigs >= {args.min_contig:,} bp are searched):\\n") print(f" {'phage':<26} {'contigs':>7} {'searched':>8} {'searched bp':>12} {'longest bp':>11} FASTA") total = 0 for name, path in genomes: recs = read_fasta(path) if not recs: print(f" {name:<26} {'—':>7} ⚠️ no sequence in {path}") continue q = genome_query(name, path, recs, args.min_contig) total += 1 print(f" {name:<26} {q['contigs']:>7} {q['searched']:>8} {q['length']:>12,} {q['longest']:>11,} {path}") extra = [] if iso: extra.append("isolation host: " + (iso.get(norm_key(name)) or "— (not in the list)")) link, flags, circ, _ = pipeline_info(name, path, args.out, args.input) circ = circ or q["circ_hdr"] extra.append("circular: " + (circ or "unknown")) if link: extra.append("coverage.html found") extra += [f"⚠️ {f}" for f in flags] if extra: print(f" {'':<26} " + "; ".join(extra)) if skipped: print(f"\\n not phage folders, skipped: {', '.join(skipped)}") if iso: names = {norm_key(n) for n, _ in genomes} unused = [k for k in iso if k not in names] print(f"\\nIsolation hosts: {len(iso)} listed" + (f"; {len(unused)} not matched to a genome folder" if unused else "")) print(f"\\nCheck only: nothing was searched. {total} genome(s) would be searched.") # ---------------------------------------------------------------- report def url_pdb(name): return "https://phagesdb.org/phages/" + urllib.parse.quote(str(name)) + "/" HTML = r""" Closest PhagesDB phages

Closest PhagesDB phages

Overview

Phages

Host Cluster Rows visible

How to read this table

Closest phage is the better match of PhagesDB and GenBank (ranked by identity × share of your genome covered, minimap2 asm20). Identity is over the aligned parts and covered is the share of your genome that aligns (like BLAST's Per. ident and Query cover, but not identical to BLASTn's numbers). Covers of closest is the share of the closest phage covered by your contigs; status: complete = one contig covers ≥90% of it, or is about as long as it (±10%) with ≥80% of yours aligned (a complete genome of a different phage), fragmented = the same with several contigs, partial = less, no relative = under 50% of your genome aligns; the grey line under each status says why. ICTV is an approximate intergenomic similarity to the closest phage (≥95% same species, ≥70% same genus). +N similar = other phages within 0.5% of the best match. All contigs of at least bp are searched. "Likely cluster" is tentative (PhagesDB assigns clusters; "?" = under 50% aligned). Click ▸ for the top hits and where your contigs align; click a bar or number in the Overview to show only those phages; click a column heading to sort. All hits and numbers: closest_phages.tsv, closest_phages_top_hits.tsv, closest_genbank_top_hits.tsv.

Figures

""" def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--input", required=True, help="folder with one subfolder (or FASTA file) per phage") ap.add_argument("--phagesdb", required=True, help="PhagesDB all-genomes FASTA (Actinobacteriophages-All.fasta)") ap.add_argument("--table", default="", help="PhagesDB data table (tab-delimited), optional") ap.add_argument("--out", required=True, help="output folder") ap.add_argument("--pattern", default="", help="FASTA file name in each subfolder, e.g. {phage}.fasta or *.fasta (default: automatic)") ap.add_argument("--threads", type=int, default=4) ap.add_argument("--top", type=int, default=5) ap.add_argument("--ncbi", choices=["none", "hosts", "all"], default="none", help="also search complete phage genomes from NCBI: by host genus, or all Caudoviricetes") ap.add_argument("--ncbi-hosts", default="Mycobacterium Gordonia Microbacterium Arthrobacter Staphylococcus", help="host genera for --ncbi hosts") ap.add_argument("--ncbi-fasta", default="", help="use this FASTA of NCBI phage genomes instead of downloading") ap.add_argument("--ncbi-dir", default="", help="where NCBI genomes are kept (default: NCBI_phages next to the PhagesDB FASTA)") ap.add_argument("--ncbi-refresh", action="store_true", help="download the NCBI genomes again") ap.add_argument("--host-genome", action="append", default=[], help="host genome FASTA (repeat for several hosts)") ap.add_argument("--checkv-db", default="", help="CheckV database folder (run CheckV when given)") ap.add_argument("--checkv-bin", default="", help="bin folder of the CheckV conda environment") ap.add_argument("--checkv-download", action="store_true", help="download the CheckV database if it is missing") ap.add_argument("--no-similarity", action="store_true", help="skip the similarity matrix") ap.add_argument("--min-contig", type=int, default=5000, help="search every contig at least this long (bp); the longest contig is used when none is") ap.add_argument("--isolation-hosts", action="append", default=[], help="phage -> isolation host list (TSV/CSV, or a PhageTracker export); repeat for several") ap.add_argument("--coverage-url", default="", help="web folder (e.g. SharePoint) that holds the phage folders with the coverage.html files; " "can also be changed in the report") ap.add_argument("--no-gc-bias", action="store_true", help="skip the read depth vs GC content analysis") ap.add_argument("--check", action="store_true", help="only list the genomes and FASTA files; no search") args = ap.parse_args() t0 = time.time() print(TOOL_VERSION) if not os.path.isdir(args.input): sys.exit(f"Input folder not found: {args.input}") if args.check: genomes, skipped = find_genomes(args.input, args.pattern, skip_dirs=[args.out]) if not genomes: sys.exit(f"No FASTA files found in the subfolders of {args.input}") check_mode(args, genomes, skipped) if not os.path.isfile(args.phagesdb): print(f"⚠️ PhagesDB FASTA not found: {args.phagesdb}") return if not os.path.isfile(args.phagesdb): sys.exit(f"PhagesDB FASTA not found: {args.phagesdb}") table = "" if args.table: for cand in (args.table, args.table + ".txt", args.table + ".tsv", args.table + ".tab"): if os.path.isfile(cand): table = cand break if not table: print(f"⚠️ PhagesDB data table not found: {args.table} (also tried .txt/.tsv/.tab); using the FASTA headers only.") os.makedirs(args.out, exist_ok=True) genomes, skipped = find_genomes(args.input, args.pattern, skip_dirs=[args.out]) if not genomes: sys.exit(f"No FASTA files found in the subfolders of {args.input}") print(f"{len(genomes)} genome(s) found in {args.input}") if NO_FASTA: print(f"⚠️ {len(NO_FASTA)} folder(s) without a genome FASTA, not included (not assembled yet, the assembly failed, " f"or it is still in the WSL work folder of a running pipeline): {', '.join(NO_FASTA)}") if skipped: print(f" not phage folders, skipped: {', '.join(skipped)}") index, meta = prepare_phagesdb(args.phagesdb, table, args.threads) by_phage = {} for pid, m in meta.items(): by_phage.setdefault(str(m.get("phage", "")).lower(), pid) # one query per genome: all its contigs >= --min-contig, joined with N spacers import tempfile qdir = tempfile.mkdtemp(prefix="closest_phages_") qpath = os.path.join(qdir, "queries.fasta") jpath = os.path.join(qdir, "genomes_joined.fasta") # for CheckV: one sequence per genome info, qoff = {}, {} with open(qpath, "w") as fh, open(jpath, "w") as fj: for i, (name, path) in enumerate(genomes): recs = read_fasta(path) if not recs: print(f"⚠️ {name}: {path} has no sequence; skipped") continue qid = f"q{i}" info[qid] = genome_query(name, path, recs, args.min_contig) write_fasta(fj, qid, info[qid].pop("seq")) byname = {r[0]: r[2] for r in recs} for k, o in enumerate(info[qid]["offs"]): write_fasta(fh, f"{qid}__c{k}", byname[o[2]]) qoff[f"{qid}__c{k}"] = (qid, o[0]) nmulti = sum(1 for q in info.values() if q["searched"] > 1) print(f"Contigs >= {args.min_contig:,} bp are searched" + (f" ({nmulti} genome(s) have more than one)" if nmulti else "")) iso = read_isolation_hosts(args.isolation_hosts) if args.isolation_hosts: print(f"Isolation hosts: {len(iso)} listed, " f"{sum(1 for q in info.values() if norm_key(q['sample']) in iso)} matched to these genomes") print(f"Searching PhagesDB with minimap2 ({args.threads} threads)…" + mem_note(), flush=True) try: hits = search(index, qpath, args.threads, qoff) except RuntimeError as e: sys.exit(str(e)) ncbi_hits, ncbi_meta, ncbi_info, ncbi_label = None, {}, {}, "" if args.ncbi != "none" or args.ncbi_fasta: hosts = [h for h in re.split(r"[,;\\s]+", args.ncbi_hosts) if h] ncbi_dir = args.ncbi_dir or os.path.join(os.path.dirname(os.path.abspath(args.phagesdb)), "NCBI_phages") try: nidx, ncbi_meta, ncbi_info = prepare_ncbi(ncbi_dir, args.ncbi, hosts, args.ncbi_fasta, args.ncbi_refresh, args.threads) ncbi_label = (os.path.basename(args.ncbi_fasta) if args.ncbi_fasta else ("all Caudoviricetes" if args.ncbi == "all" else "phages of " + ", ".join(hosts))) print(f"Searching {len(ncbi_meta):,} NCBI genomes ({ncbi_label}" f"{'; downloaded ' + ncbi_info['downloaded'] if ncbi_info.get('downloaded') else ''})…" + mem_note(), flush=True) ncbi_hits = search(nidx, qpath, args.threads, qoff) except Exception as e: # keep the PhagesDB results even if NCBI fails print(f"⚠️ NCBI search skipped: {e}") ncbi_hits = None # ---- genome checks hostres, checkv, sim, pair, checks = {}, {}, {}, {}, {"host": False, "checkv": False, "sim": False, "host_names": []} hosts_ok = [h for h in args.host_genome if h and os.path.isfile(h)] for h in args.host_genome: if h and not os.path.isfile(h): print(f"⚠️ Host genome not found: {h}") if hosts_ok: print(f"Checking for host DNA ({len(hosts_ok)} host genome(s))…" + mem_note(), flush=True) try: hostres = host_check(qpath, hosts_ok, args.threads, qoff, {k: v["length"] for k, v in info.items()}) checks["host"] = True checks["host_names"] = [re.sub(r"\\.(fasta|fa|fna|fas)(\\.gz)?$", "", os.path.basename(h), flags=re.I) for h in hosts_ok] except RuntimeError as e: print(f"⚠️ Host check skipped: {e}") if args.checkv_db: try: checkv = run_checkv(jpath, args.out, args.checkv_bin, args.checkv_db, args.checkv_download, args.threads) checks["checkv"] = True except Exception as e: print(f"⚠️ CheckV skipped: {e}") if not args.no_similarity and len(info) > 1: print("Comparing all genomes with each other (intergenomic similarity)…" + mem_note(), flush=True) try: sim, pair = similarity_matrix(qpath, info, args.threads, qoff) checks["sim"] = True except RuntimeError as e: print(f"⚠️ Similarity matrix skipped: {e}") species = {q: g for g in groups_at(list(info), sim, 95.0) for q in g} if checks["sim"] else {} genus = {q: g for g in groups_at(list(info), sim, 70.0) for q in g} if checks["sim"] else {} gnum, snum = {}, {} for q in sorted(info, key=lambda x: info[x]["sample"].lower()): if genus: gnum.setdefault(tuple(genus[q]), len(gnum) + 1) snum.setdefault(tuple(species[q]), len(snum) + 1) import shutil shutil.rmtree(qdir, ignore_errors=True) # phages in this run that are the same phage: >= 98% similarity, or >= 98% identity over >= 80% of the shorter # genome (a fragmented assembly can't reach 98% similarity); the value shown is the identity dups = {} for (a, b), v in sim.items(): idn, cs = pair.get((a, b), (0.0, 0.0)) if v >= 98 or (idn >= 98 and cs >= 80): dups.setdefault(a, []).append([info[b]["sample"], idn if v < 98 else v]) dups.setdefault(b, []).append([info[a]["sample"], idn if v < 98 else v]) if not args.no_gc_bias: print("Read depth vs GC content (from the pipeline's coverage tables)…", flush=True) rows, long_rows, ncbi_long = [], [], [] clu_len, name_clu = {}, {} for _tid, _m in meta.items(): try: _ln = int(str(_m.get("genome_length_table") or _m.get("length") or "0").replace(",", "")) except ValueError: _ln = 0 if _m.get("cluster") and _ln > 1000: clu_len.setdefault(_m["cluster"], []).append(_ln) if _m.get("phage"): name_clu[_m["phage"].lower()] = _m.get("cluster", "") for qid, q in info.items(): L = q["length"] pr = rank_hits(hits.get(qid), L) nr = rank_hits(ncbi_hits.get(qid), L) if ncbi_hits is not None else [] row = {"sample": q["sample"], "file": q["file"], "contigs": q["contigs"], "searched": q["searched"], "longest": q["longest"], "contig_list": [[o[2], o[3]] for o in q["offs"]], "length": L, "gc": q["gc"], "hit": "", "hit_url": "", "self": False, "cluster": "", "identity": "", "aligned": "", "len_diff": "", "host": "", "temperate": "", "morphotype": "", "ends": "", "accession": "", "gc_hit": "", "other": "", "other_url": "", "other_cluster": "", "other_identity": "", "other_aligned": "", "next": [], "note": "", "ncbi": None, "ncbi_next": [], "host_pct": "", "host_name": "", "cv_quality": "", "cv_completeness": "", "cv_method": "", "cv_contamination": "", "cv_viral_genes": "", "cv_host_genes": "", "cv_provirus": "", "cv_warnings": "", "species": "", "species_with": [], "genus": "", "genus_with": [], "best": {}, "status": "", "ties": [], "dups": sorted(dups.get(qid, []), key=lambda x: -x[1]), "iso_host": iso.get(norm_key(q["sample"]), ""), "host_mismatch": "", "cov_html": "", "pipe_flags": [], "top_pdb": [], "top_gb": []} row["cov_html"], row["pipe_flags"], row["circular"], row["cov_rel"] = pipeline_info(q["sample"], q["file"], args.out, args.input) row["circular"] = row["circular"] or q["circ_hdr"] row.update(pipeline_reads(q["file"])) sq = sequenced(os.path.dirname(os.path.abspath(q["file"]))) row["seq_reads"], row["seq_bases"], row["seq_src"] = sq["reads"], sq["bases"], sq["src"] row["gcb"] = gc_bias(q["file"], q["sample"]) if not args.no_gc_bias else None row["contig_hits"] = contig_hits(hits.get(qid), q["offs"], meta) if checks["host"]: hp = hostres.get(qid) row["host_pct"], row["host_name"] = (hp[0], hp[1]) if hp else (0.0, "") if checks["checkv"] and qid in checkv: c = checkv[qid] row.update({"cv_quality": c.get("checkv_quality", ""), "cv_completeness": c.get("completeness", "NA"), "cv_method": c.get("completeness_method", ""), "cv_contamination": c.get("contamination", ""), "cv_viral_genes": c.get("viral_genes", ""), "cv_host_genes": c.get("host_genes", ""), "cv_provirus": c.get("provirus", ""), "cv_warnings": c.get("warnings", "")}) if checks["sim"]: row["species"] = f"S{snum[tuple(species[qid])]}" row["genus"] = f"G{gnum[tuple(genus[qid])]}" row["species_with"] = [info[x]["sample"] for x in species[qid] if x != qid] row["genus_with"] = [info[x]["sample"] for x in genus[qid] if x != qid] if ncbi_hits is not None: for rank, x in enumerate(nr[:args.top], 1): m = ncbi_meta.get(x["id"], {}) ncbi_long.append([q["sample"], rank, x["id"], m.get("description", ""), m.get("host", ""), x["ident"], x["qcov"], x["aln"], x["tlen"], x["scov"], L - x["tlen"], url_ncbi(x["id"])]) row["top_gb"].append([x["id"], m.get("description", ""), m.get("host", ""), x["ident"], x["qcov"], x["scov"], x["tlen"], url_ncbi(x["id"])]) if nr: x = nr[0] m = ncbi_meta.get(x["id"], {}) row["ncbi"] = {"acc": x["id"], "desc": m.get("description", ""), "host": m.get("host", ""), "ident": x["ident"], "qcov": x["qcov"], "aln": x["aln"], "slen": x["tlen"], "scov": x["scov"], "len_diff": f"{L - x['tlen']:+d}", "url": url_ncbi(x["id"])} row["ncbi_next"] = [[y["id"], ncbi_meta.get(y["id"], {}).get("description", ""), y["ident"], y["qcov"], url_ncbi(y["id"])] for y in nr[1:4]] else: row["ncbi"] = {} for rank, x in enumerate(pr[:args.top], 1): m = meta.get(x["id"], {}) nm = m.get("phage") or x["id"] long_rows.append([q["sample"], rank, nm, m.get("cluster", ""), x["ident"], x["qcov"], x["tlen"], L - x["tlen"], m.get("host", ""), m.get("temperate", ""), m.get("morphotype", ""), m.get("ends", ""), m.get("accession", ""), url_pdb(nm)]) row["top_pdb"].append([nm, m.get("cluster", ""), m.get("host", ""), x["ident"], x["qcov"], x["scov"], x["tlen"], url_pdb(nm)]) if pr: x = pr[0] m = meta.get(x["id"], {}) name = m.get("phage") or x["id"] row.update({"hit": name, "hit_url": url_pdb(name), "self": name.lower() == q["sample"].lower(), "cluster": m.get("cluster", ""), "identity": x["ident"], "aligned": x["qcov"], "len_diff": f"{L - x['tlen']:+d}", "host": m.get("host", ""), "temperate": m.get("temperate", ""), "morphotype": m.get("morphotype", ""), "ends": m.get("ends", ""), "accession": m.get("accession", ""), "gc_hit": m.get("gc_pct", "")}) for y in pr: n2 = meta.get(y["id"], {}).get("phage") or y["id"] if n2.lower() != q["sample"].lower(): row.update({"other": n2, "other_url": url_pdb(n2), "other_cluster": meta.get(y["id"], {}).get("cluster", ""), "other_identity": y["ident"], "other_aligned": y["qcov"]}) break row["next"] = [[meta.get(y["id"], {}).get("phage") or y["id"], meta.get(y["id"], {}).get("cluster", ""), y["ident"], y["qcov"], url_pdb(meta.get(y["id"], {}).get("phage") or y["id"])] for y in pr[1:4]] if x["qcov"] < 50: row["note"] = "no close relative (<50% aligned)" elif row["self"]: row["note"] = "same name as the PhagesDB phage" else: row["note"] = "no PhagesDB phage aligns" # closest phage overall: the better of the PhagesDB and GenBank top hits (identity x query cover) cands = [] if pr: x = pr[0] cands.append({"score": float(x["ident"]) * float(x["qcov"]), "db": "PhagesDB", "name": row["hit"], "url": row["hit_url"], "extra": row["cluster"], "host": row["host"], "acc": row["accession"], "x": x, "d": hits[qid][x["id"]]}) if nr: x, nb = nr[0], row["ncbi"] short = re.sub(r",?\\s*(complete|partial)\\s+(genome|sequence|cds).*$", "", nb["desc"] or nb["acc"], flags=re.I).strip() cands.append({"score": float(x["ident"]) * float(x["qcov"]), "db": "GenBank", "name": short or nb["acc"], "url": nb["url"], "extra": nb["acc"], "host": nb["host"], "acc": nb["acc"], "x": x, "d": ncbi_hits[qid][x["id"]]}) if cands: b = max(cands, key=lambda c: c["score"]) x = b["x"] sv = ictv_sim(x["ident"], x["aln"], x["taln"], L, x["tlen"]) row["best"] = {"name": b["name"], "url": b["url"], "db": b["db"], "ident": x["ident"], "qcov": x["qcov"], "extra": b["extra"], "len_diff": f"{L - x['tlen']:+d}", "host": b["host"], "acc": b["acc"], "slen": x["tlen"], "scov": x["scov"], "sim": sv, "rank": ictv_rank(sv), "rank_note": ("" if float(x["scov"]) >= 90 or L >= 0.9 * x["tlen"] else "your contigs cover under 90% of the closest phage, so the similarity is underestimated"), "blocks": aln_blocks(b["d"], q["offs"])} if b["db"] == "GenBank" and row["note"].startswith("no "): row["note"] = "closest relative is in GenBank, not PhagesDB" # near-equal matches: identity x cover within 0.5% of the best, in either database thr, seen, ties = b["score"] * 0.995, {b["name"].lower(), str(b["acc"]).split(".")[0].lower()}, [] for db, lst in (("PhagesDB", pr), ("GenBank", nr)): for y in lst: sc = float(y["ident"]) * float(y["qcov"]) if sc < thr: break if db == "PhagesDB": m = meta.get(y["id"], {}) nm, acc, url = m.get("phage") or y["id"], m.get("accession", ""), url_pdb(m.get("phage") or y["id"]) else: m = ncbi_meta.get(y["id"], {}) nm = re.sub(r",?\\s*(complete|partial)\\s+(genome|sequence|cds).*$", "", m.get("description", "") or y["id"], flags=re.I).strip() or y["id"] acc, url = y["id"], url_ncbi(y["id"]) keys = {nm.lower(), str(acc).split(".")[0].lower()} - {""} if keys & seen: continue seen |= keys ties.append([nm, db, y["ident"], y["qcov"], url]) row["ties"] = ties[:20] n_ties = len(ties) else: row["best"] = {} n_ties = 0 bst = row["best"] row["status_note"] = "" # v3.1: judge the genome on its own contigs (those that match the closest phage), not on all contigs together offs = q["offs"] chosen, Lm, sm, qm = list(range(len(offs))), L, 0.0, 0.0 if bst: chosen, sm, qm, Lm, qa, ta = main_contigs(b["d"], offs, x["tlen"]) if not chosen: chosen, Lm = list(range(len(offs))), L bst["qcov_all"], bst["scov_all"] = bst["qcov"], bst["scov"] bst["qcov"], bst["scov"] = qm, sm bst["len_diff"] = f"{Lm - x['tlen']:+d}" if chosen and qa: sv2 = ictv_sim(x["ident"], qa, ta, Lm, x["tlen"]) bst["sim"], bst["rank"] = sv2, ictv_rank(sv2) bst["rank_note"] = ("" if sm >= 90 or Lm >= 0.9 * x["tlen"] else "your contigs cover under 90% of the closest phage, so the similarity is underestimated") bst["blocks"] = [bl for bl in bst["blocks"] if bl[2] in chosen] others = [k for k in range(len(offs)) if k not in chosen] row["main_contigs"] = [offs[k][2] for k in chosen] row["main_len"] = Lm row["other_contigs"] = [len(others), sum(offs[k][3] for k in others)] cc = cluster_call(hits.get(qid), offs, chosen, Lm, meta, clu_len) row["cluster_call"] = cc row["cluster_letter"] = cc["letter"] if cc["hit"]: row["hit"], row["hit_url"], row["cluster"], row["identity"], row["aligned"] = (cc["hit"], url_pdb(cc["hit"]), cc["cluster"], cc["ident"], cc["span"]) elif row["hit"]: row["cluster"], row["aligned"] = "", 0 row["gb_group"] = gb_group(bst.get("name", "")) if bst and bst.get("db") == "GenBank" and not cc["letter"] else "" if bst and bst.get("db") == "GenBank": nm = re.sub(r"^.*?\\b(phage|virus)\\s+", "", bst["name"]).strip().lower() bst["pdb_cluster"] = name_clu.get(nm, "") n_main = len(chosen) rng = cc["range"] if not bst or qm < 50: row["status"] = "no relative" elif sm >= 90: row["status"] = "complete" if n_main == 1 else "fragmented" elif 0.9 <= (Lm / float(x["tlen"])) <= 1.1 and qm >= 80: # v2.2: about as long as the closest phage and mostly aligned, but only part of the relative matches: # a complete genome of a different phage, not a missing piece row["status"] = "complete" if n_main == 1 else "fragmented" row["status_note"] = (f"{Lm:,} bp, about as long as the closest phage ({int(x['tlen']):,} bp), but only " f"{sm}% of it matches: likely complete, and a different phage") elif n_main == 1 and rng and 0.97 * rng[0] <= Lm <= 1.03 * rng[1] and cc["span"] >= 50: # v3.1: one contig as long as the phages of its cluster: a complete genome of a new member of the cluster row["status"] = "complete" row["status_note"] = (f"{Lm:,} bp, within cluster {cc['letter']}'s range ({rng[0] / 1000:.1f}–{rng[1] / 1000:.1f} kb), " f"{cc['ident']}% identical to its closest relatives: a complete genome of a new {cc['cluster'].rstrip('?')} phage") else: row["status"] = "partial" cv, ev, sc = completeness(n_main, row["circular"], Lm, rng, sm if bst else "", row["cv_quality"], cc["letter"] or "?") row["completeness"], row["completeness_ev"], row["completeness_score"] = cv, ev, sc row["contains"], row["pair"] = [], {} ig, bg = host_genus(row["iso_host"]), host_genus(bst.get("host", "")) if ig and bg and ig.lower() != bg.lower() and row["status"] != "no relative": row["host_mismatch"] = f"isolated on {ig}, closest phage infects {bg}" row["_ties_n"] = n_ties rows.append(row) post_rows(rows, checks) rows.sort(key=lambda x: x["sample"].lower()) cols = ["phage", "genome_file", "contigs", "length_bp", "gc_pct", "closest_phage", "closest_source", "closest_identity_pct", "closest_query_cover_pct", "closest_link", "closest_phagesdb", "likely_cluster", "identity_pct", "aligned_pct", "length_difference_bp", "host", "temperate", "morphotype", "ends", "accession", "phagesdb_gc_pct", "closest_other_phage", "other_cluster", "other_identity_pct", "other_aligned_pct", "note", "phagesdb_link", "genbank_accession", "genbank_description", "genbank_host", "genbank_per_ident", "genbank_query_cover", "genbank_aligned_bp", "genbank_subject_length", "genbank_subject_cover", "genbank_length_difference_bp", "genbank_link", "host_dna_pct", "host_matched", "checkv_quality", "checkv_completeness_pct", "checkv_completeness_method", "checkv_contamination_pct", "checkv_viral_genes", "checkv_host_genes", "checkv_provirus", "checkv_warnings", "species_group_95pct", "same_species_as", "genus_group_70pct", "same_genus_as", "contigs_searched", "closest_length_bp", "closest_covered_pct", "status", "similar_matches", "ictv_similarity_to_closest_pct", "ictv_rank", "same_phage_in_run", "isolation_host", "host_mismatch", "coverage_html", "pipeline_flags", "circular", "status_note", "bases_sequenced", "median_read_depth", "pct_bases_low_depth", "low_depth_regions", "gc_bias_rho", "gc_bias_log2_depth_per_10pct_gc", "depth_in_windows_5pct_below_mean_gc", "genome_contigs", "genome_length_bp", "other_contigs", "other_contigs_bp", "cluster_agreement", "cluster_note", "completeness", "completeness_evidence", "contains_dna_of", "ligase_pair", "closest_query_cover_all_contigs_pct"] with open(os.path.join(args.out, "closest_phages.tsv"), "w") as out: out.write("\\t".join(cols) + "\\n") for x in rows: out.write("\\t".join(str(v) for v in [ x["sample"], x["file"], x["contigs"], x["length"], x["gc"]] + [x["best"].get(k, "") for k in ("name", "db", "ident", "qcov", "url")] + [x["hit"], x["cluster"], x["identity"], x["aligned"], x["len_diff"], x["host"], x["temperate"], x["morphotype"], x["ends"], x["accession"], x["gc_hit"], x["other"], x["other_cluster"], x["other_identity"], x["other_aligned"], x["note"], x["hit_url"]] + [(x["ncbi"] or {}).get(k, "") for k in ("acc", "desc", "host", "ident", "qcov", "aln", "slen", "scov", "len_diff", "url")] + [x["host_pct"], x["host_name"], x["cv_quality"], x["cv_completeness"], x["cv_method"], x["cv_contamination"], x["cv_viral_genes"], x["cv_host_genes"], x["cv_provirus"], x["cv_warnings"], x["species"], ", ".join(x["species_with"]), x["genus"], ", ".join(x["genus_with"])] + [x["searched"], x["best"].get("slen", ""), x["best"].get("scov", ""), x["status"], ", ".join(t[0] for t in x["ties"]), x["best"].get("sim", ""), x["best"].get("rank", "") if x["best"] else "", ", ".join(f"{n} ({v}%)" for n, v in x["dups"]), x["iso_host"], x["host_mismatch"], x["cov_html"], "; ".join(x["pipe_flags"]), x["circular"], x["status_note"], x["bases"], x["med_depth"], x["pct_low"], x["dips"]] + ([x["gcb"]["rho"], x["gcb"]["slope"], x["gcb"]["low5"] if x["gcb"]["low5"] is not None else ""] if x.get("gcb") else ["", "", ""]) + [", ".join(x["main_contigs"]), x["main_len"], x["other_contigs"][0], x["other_contigs"][1], x["cluster_call"]["agree"], x["cluster_call"]["note"], x["completeness"], "; ".join(("yes: " if e[1] else ("no: " if e[1] is False else "unknown: ")) + e[0] for e in x["completeness_ev"]), ", ".join(f"{n} ({lab})" for n, v, lab in x["contains"]), (("better than " if x["pair"]["better"] else "worse than ") + x["pair"]["with"] + " (" + x["pair"]["why"] + ")") if x["pair"] else "", x["best"].get("qcov_all", "") if x["best"] else ""]) + "\\n") with open(os.path.join(args.out, "closest_phages_top_hits.tsv"), "w") as out: out.write("phage\\trank\\tphagesdb_phage\\tcluster\\tidentity_pct\\taligned_pct\\tphagesdb_length_bp\\tlength_difference_bp" "\\thost\\ttemperate\\tmorphotype\\tends\\taccession\\tphagesdb_link\\n") for lr in long_rows: out.write("\\t".join(str(v) for v in lr) + "\\n") if ncbi_hits is not None: with open(os.path.join(args.out, "closest_genbank_top_hits.tsv"), "w") as out: out.write("phage\\trank\\taccession\\tdescription\\thost\\tper_ident\\tquery_cover\\taligned_bp\\tsubject_length" "\\tsubject_cover\\tlength_difference_bp\\tgenbank_link\\n") for lr in ncbi_long: out.write("\\t".join(str(v) for v in lr) + "\\n") matrix = None if checks["sim"]: order = sorted(info, key=lambda q: (gnum[tuple(genus[q])], snum[tuple(species[q])], info[q]["sample"].lower())) names = [info[q]["sample"] for q in order] vals = [[100.0 if a == b else sim.get(tuple(sorted((a, b))), 0.0) for b in order] for a in order] with open(os.path.join(args.out, "similarity_matrix.tsv"), "w") as out: out.write("phage\\t" + "\\t".join(names) + "\\n") for n, row_v in zip(names, vals): out.write(n + "\\t" + "\\t".join(str(v) for v in row_v) + "\\n") matrix = {"names": names, "vals": vals, "species": [f"S{snum[tuple(species[q])]}" for q in order], "genus": [f"G{gnum[tuple(genus[q])]}" for q in order]} data = {"rows": rows, "input": args.input, "pdb_file": os.path.basename(args.phagesdb), "pdb_n": len(meta), "checks": checks, "matrix": matrix, "when": time.strftime("%Y-%m-%d %H:%M"), "ncbi": ncbi_hits is not None, "ncbi_n": len(ncbi_meta), "ncbi_label": ncbi_label, "ncbi_date": ncbi_info.get("downloaded", ""), "ncbi_max_age": NCBI_MAX_AGE_DAYS, "min_contig": args.min_contig, "cov_base": args.coverage_url, "has_iso": bool(iso), "version": TOOL_VERSION, "no_fasta": [dict(sample=n, **sequenced(os.path.join(args.input, n))) for n in NO_FASTA]} page = HTML.replace("__DATA__", json.dumps(data, separators=(",", ":")).replace("= 50) print(f"\\n{len(rows)} genomes: {n_close} with a close PhagesDB relative, {len(rows) - n_close} without " f"({time.time() - t0:.0f} s)") if checks["host"]: nh = sum(1 for x in rows if x["host_pct"] != "" and float(x["host_pct"]) >= 50) print(f"Host DNA: {nh} genome(s) are mostly host DNA (≥50% matches a host genome)") if checks["checkv"]: qual = {} for x in rows: qual[x["cv_quality"] or "not run"] = qual.get(x["cv_quality"] or "not run", 0) + 1 print("CheckV: " + ", ".join(f"{k} {v}" for k, v in sorted(qual.items()))) if checks["sim"]: multi = [x for x in rows if x["species_with"]] print(f"Similarity: {len(snum)} species group(s) (≥95%), {len(gnum)} genus group(s) (≥70%); " f"{len(multi)} genome(s) share a species group with another") print(f"Matrix: {os.path.join(args.out, 'similarity_matrix.tsv')}") st = {} for x in rows: st[x["status"]] = st.get(x["status"], 0) + 1 print("Status: " + ", ".join(f"{k} {v}" for k, v in sorted(st.items()))) gb = [x for x in rows if x.get("gcb")] if gb: neg = sum(1 for x in gb if x["gcb"]["rho"] >= 0.3) print(f"GC bias: {len(gb)} phage(s) with coverage tables; {neg} with clearly lower depth in AT-rich DNA (rho >= 0.3), " f"{sum(1 for x in gb if x['gcb']['rho'] <= -0.3)} with lower depth in GC-rich DNA (rho <= -0.3)") nd = sum(1 for x in rows if x["dups"]) if nd: print(f"Same phage in this run (>= 98% identity over >= 80% of the shorter genome, or >= 98% similar): {nd} genome(s)") print(f"Table: {os.path.join(args.out, 'closest_phages.tsv')}") print(f"Report: {os.path.join(args.out, 'closest_phages.html')}") if __name__ == "__main__": main() `; // v340: shared with the Nanopore phage pipeline cell on the same page (it uses this copy when it is newer than its own) const CLOSEST_PY_VERSION = 356; try { if (!window.__closestPhagesPy || window.__closestPhagesPy.version <= CLOSEST_PY_VERSION) window.__closestPhagesPy = { version: CLOSEST_PY_VERSION, code: CLOSEST_PY, from: "Closest phages generator" }; } catch (_) {} function readSettings() { const clean = (v, d) => (String(v || "").trim() || d).replace(/[\t\r\n]+/g, " "); return { INPUT_DIR: toPath($("inputDir").value), PATTERN: clean($("pattern").value, ""), TOP: Math.min(20, Math.max(1, Math.floor(Number($("top").value) || 5))), MIN_CONTIG: Math.max(500, Math.floor(Number($("minContig").value) || 5000)), ISO_HOSTS_FILE: toPath($("isoHostsFile").value), PDB_FASTA: toPath($("pdbFasta").value), PDB_TABLE: toPath($("pdbTable").value), OUT_DIR: toPath($("outDir").value), CONDA_ENV: clean($("condaEnv").value, ""), CONDA_BASE_INPUT: String($("condaBase").value || "").trim(), THREADS: Math.max(1, Math.floor(Number($("threads").value) || 8)), NCBI_MODE: $("ncbiMode").value, NCBI_HOSTS: [...new Set([...[...wrap.querySelectorAll(".cp-genus:checked")].map(x => x.value), ...clean($("ncbiHostsOther").value, "").split(/[\s,;]+/).filter(Boolean)])].join(" "), NCBI_FASTA: toPath($("ncbiFasta").value), NCBI_REFRESH: $("ncbiRefresh").checked ? "true" : "false", HOST_GENOME_1: toPath($("hostGenome1").value), HOST_GENOME_2: toPath($("hostGenome2").value), HOST_GENOME_3: toPath($("hostGenome3").value), HOST_GENOME_4: toPath($("hostGenome4").value), HOST_GENOME_5: toPath($("hostGenome5").value), RUN_CHECKV: $("runCheckv").checked ? "true" : "false", CHECKV_ENV: clean($("checkvEnv").value, "checkv"), CHECKV_INSTALL: $("checkvInstall").checked ? "true" : "false", CHECKV_DB: clean($("checkvDb").value, "~/databases/checkv"), CHECKV_DOWNLOAD: $("checkvDownload").checked ? "true" : "false", RUN_SIMILARITY: $("runSimilarity").checked ? "true" : "false", OPEN_REPORT: $("openReport").checked ? "true" : "false" }; } function validate(s) { const errors = [], warnings = []; if (!s.INPUT_DIR) errors.push("Enter the folder with the phage genomes."); if (/\.(fastq|fq|pod5|fasta|fa|fna|html|tsv)(\.gz)?$/i.test(s.INPUT_DIR)) errors.push("Genome folder: enter the folder, not a file."); if (!s.PDB_FASTA) errors.push("Enter the PhagesDB genomes FASTA."); else if (!/\.(fa|fasta|fna|fas)(\.gz)?$/i.test(s.PDB_FASTA)) warnings.push("The PhagesDB FASTA path does not end in .fasta, .fa or .fna; check the file name."); for (const [k, label] of [["INPUT_DIR", "Genome folder"], ["PDB_FASTA", "PhagesDB FASTA"], ["PDB_TABLE", "Data table"], ["OUT_DIR", "Report folder"]]) { if (/^\\\\/.test(s[k])) errors.push(`${label}: network (UNC) paths can't be converted; map the drive to a letter first.`); } if (/(^|\/)phage_fastq(\/|$)/i.test(s.INPUT_DIR)) warnings.push("That looks like the FASTQ folder; use a folder with assembled genomes (e.g. the pipeline's output folder)."); if (s.NCBI_MODE === "hosts" && !s.NCBI_HOSTS && !s.NCBI_FASTA) errors.push("Tick at least one host genus for the GenBank search (or choose No)."); if (s.NCBI_MODE === "hosts" && /[^A-Za-z .-]/.test(s.NCBI_HOSTS)) warnings.push("Other genera: use genus names such as Vibrio Clostridium."); if (s.NCBI_MODE === "all" && !s.NCBI_FASTA) warnings.push("All Caudoviricetes is a large download (several GB) and its index needs a lot of memory; host genera is usually enough."); if (s.RUN_CHECKV === "true" && !/^[A-Za-z0-9._-]+$/.test(s.CHECKV_ENV)) errors.push("The CheckV environment name contains unsupported characters."); if (!isMac() && s.RUN_CHECKV === "true" && /^[A-Za-z]:/.test($("checkvDb").value.trim())) warnings.push("CheckV database on a Windows drive works but is slower; a WSL folder such as ~/databases/checkv is faster."); for (const i of [1, 2, 3, 4, 5]) { const g = s[`HOST_GENOME_${i}`]; if (g && !/\.(fa|fasta|fna|fas)(\.gz)?$/i.test(g)) warnings.push(`Host genome ${i}: the path does not end in .fasta, .fa or .fna.`); } if (s.ISO_HOSTS_FILE && !/\.(csv|tsv|txt|tab)$/i.test(s.ISO_HOSTS_FILE)) warnings.push("Isolation host file: use a .csv, .tsv or .txt file."); if (isoText() && !isoText().split("\n").every((l) => /[\t,;:]|\s{2,}/.test(l.trim()) || /^\S+\s+[A-Z]/.test(l.trim()))) warnings.push("Isolation hosts: put the phage and the host on each line, separated by a tab or comma."); if (isMac()) { const winP = ["inputDir", "pdbFasta", "pdbTable", "outDir", "ncbiFasta", "isoHostsFile"].filter((id) => $(id) && /^[A-Za-z]:[\\/]|^\/mnt\/[a-z]\//.test($(id).value.trim())); if (winP.length) warnings.push("These look like Windows paths, but this computer is a Mac: use Mac folders such as ~/Nanopore/Results (" + winP.map((id) => ({ inputDir: "genome folder", pdbFasta: "PhagesDB FASTA", pdbTable: "PhagesDB table", outDir: "report folder", ncbiFasta: "GenBank FASTA", isoHostsFile: "isolation host file" })[id]).join(", ") + ")."); } if (!s.CONDA_ENV) warnings.push("No conda environment: minimap2 must already be on the PATH" + (isMac() ? "." : " in WSL.")); else if (!/^(\/[^'\s]+|[A-Za-z0-9._-]+)$/.test(s.CONDA_ENV)) errors.push("The conda environment name contains unsupported characters."); return { errors, warnings }; } // pasted isolation hosts: phagehost lines (tabs, commas, semicolons or colons kept as typed) function isoText() { return String($("isoHosts").value || "").split(/\r?\n/).map((l) => l.replace(/\s+$/, "")).filter((l) => l.trim() && l.trim() !== "ISO_HOSTS_END").join("\n"); } function makeScript(s) { const settings = Object.entries(s).map(([k, v]) => `${k}=${typeof v === "number" ? v : shQuote(v)}`).join("\n"); const iso = isoText(); return [ "#!/usr/bin/env bash", "# Closest PhagesDB and GenBank phages for a folder of genomes (generated in Observable)", "# Run with: bash closest_phages.sh", "# Generator: closest-phages v356 (2026-11-06) for " + (isMac() ? "macOS" : "Windows (WSL)"), "# Check first (lists the phages and the FASTA chosen, no search): bash closest_phages.sh --check", "", "(", "set -Eeuo pipefail", "", "# ------------------------- settings -------------------------", settings, "# -------------------------------------------------------------", 'CHECK_ONLY=false; [[ "${1:-}" == "--check" ]] && CHECK_ONLY=true', '# ~ in a folder path means the home folder', 'for _v in INPUT_DIR PDB_FASTA PDB_TABLE OUT_DIR NCBI_FASTA ISO_HOSTS_FILE HOST_GENOME_1 HOST_GENOME_2 HOST_GENOME_3 HOST_GENOME_4 HOST_GENOME_5; do _x="${!_v}"; if [[ "$_x" == "~/"* ]]; then printf -v "$_v" "%s" "$HOME/${_x:2}"; fi; done', "", 'fail() { echo "❌ $*" >&2; exit 1; }', '[[ -d "$INPUT_DIR" ]] || fail "Genome folder not found: $INPUT_DIR"', '[[ -s "$PDB_FASTA" ]] || fail "PhagesDB FASTA not found: $PDB_FASTA"', 'if [[ -z "$OUT_DIR" ]]; then OUT_DIR="$INPUT_DIR/closest_phages_$(date +%Y%m%d_%H%M%S)"; fi', 'mkdir -p "$OUT_DIR"', '# progress is also written to a log in the report folder (kept if the terminal stops)', '[[ "$CHECK_ONLY" == "true" ]] || { exec > >(tee -a "$OUT_DIR/closest_phages_run.log") 2>&1; echo "Log: $OUT_DIR/closest_phages_run.log"; }', "", "# ---- conda environment with minimap2 ----", 'if [[ -n "$CONDA_ENV" ]]; then', ' CONDA_BASE=""', ' for c in "${CONDA_BASE_INPUT/#\\~/$HOME}" "$(conda info --base 2>/dev/null || true)" "$HOME/miniconda" "$HOME/miniconda3" "$HOME/anaconda3" "$HOME/miniforge3" "$HOME/mambaforge"; do', ' if [[ -n "$c" && -f "$c/etc/profile.d/conda.sh" ]]; then CONDA_BASE="$c"; break; fi', " done", ' [[ -n "$CONDA_BASE" ]] || fail "Conda was not found. Enter the conda folder in the generator (e.g. ~/miniconda)."', ' # shellcheck disable=SC1091', ' source "$CONDA_BASE/etc/profile.d/conda.sh"', ' set +u; conda activate "$CONDA_ENV" || { set -u; fail "Could not activate the conda environment: $CONDA_ENV"; }; set -u', "fi", 'command -v minimap2 >/dev/null 2>&1 || fail "minimap2 is not available. Use the Medaka environment (it has minimap2) or install it: conda install -c bioconda minimap2"', 'echo "minimap2: $(command -v minimap2) ($(minimap2 --version))"', "", "# ---- CheckV (its own conda environment) ----", 'CHECKV_ARGS=()', 'if [[ "$RUN_CHECKV" == "true" && "$CHECK_ONLY" != "true" ]]; then', ' CHECKV_DB_DIR="${CHECKV_DB/#\\~/$HOME}"', ' CHECKV_BIN=""', ' if [[ -n "${CONDA_BASE:-}" && -x "$CONDA_BASE/envs/$CHECKV_ENV/bin/checkv" ]]; then', ' CHECKV_BIN="$CONDA_BASE/envs/$CHECKV_ENV/bin"', ' elif command -v checkv >/dev/null 2>&1; then', ' CHECKV_BIN="$(dirname "$(command -v checkv)")"', ' elif [[ "$CHECKV_INSTALL" == "true" && -n "${CONDA_BASE:-}" ]]; then', ' echo "Creating the conda environment \"$CHECKV_ENV\" with CheckV (once; a few minutes)…"', ' if command -v mamba >/dev/null 2>&1; then PM=mamba; else PM=conda; fi', ' if "$PM" create -y -n "$CHECKV_ENV" -c conda-forge -c bioconda checkv && [[ -x "$CONDA_BASE/envs/$CHECKV_ENV/bin/checkv" ]]; then CHECKV_BIN="$CONDA_BASE/envs/$CHECKV_ENV/bin"; fi', ' fi', ' if [[ -n "$CHECKV_BIN" ]]; then', ' CHECKV_ARGS=(--checkv-db "$CHECKV_DB_DIR" --checkv-bin "$CHECKV_BIN")', ' [[ "$CHECKV_DOWNLOAD" == "true" ]] && CHECKV_ARGS+=(--checkv-download)', ' echo "CheckV: $CHECKV_BIN/checkv (database: $CHECKV_DB_DIR)"', ' else', ' echo "⚠️ CheckV is not installed (environment \"$CHECKV_ENV\"); continuing without it. Install: conda create -n $CHECKV_ENV -c conda-forge -c bioconda checkv"', ' fi', 'fi', 'HOST_ARGS=()', 'for hg in "$HOST_GENOME_1" "$HOST_GENOME_2" "$HOST_GENOME_3" "$HOST_GENOME_4" "$HOST_GENOME_5"; do [[ -n "$hg" ]] && HOST_ARGS+=(--host-genome "$hg"); done', 'SIM_ARGS=(); [[ "$RUN_SIMILARITY" == "true" ]] || SIM_ARGS=(--no-similarity)', 'ISO_ARGS=()', 'if [[ -n "$ISO_HOSTS_FILE" ]]; then [[ -s "$ISO_HOSTS_FILE" ]] && ISO_ARGS+=(--isolation-hosts "$ISO_HOSTS_FILE") || echo "⚠️ Isolation host file not found: $ISO_HOSTS_FILE"; fi', ...(iso ? ['cat > "$OUT_DIR/isolation_hosts_pasted.txt" <<\'ISO_HOSTS_END\'', iso, "ISO_HOSTS_END", 'ISO_ARGS+=(--isolation-hosts "$OUT_DIR/isolation_hosts_pasted.txt")'] : []), 'CHECK_ARGS=(); [[ "$CHECK_ONLY" == "true" ]] && CHECK_ARGS=(--check)', "", "# ---- the search tool (embedded) ----", 'TOOL_DIR="$(mktemp -d)"', 'trap \'rm -rf "$TOOL_DIR"\' EXIT', 'cat > "$TOOL_DIR/closest_phages.py" <<\'CLOSEST_PHAGES_PY\'', CLOSEST_PY, "CLOSEST_PHAGES_PY", "", 'echo "Genomes: $INPUT_DIR"', 'echo "PhagesDB: $PDB_FASTA"', 'echo "Report: $OUT_DIR"', 'python3 "$TOOL_DIR/closest_phages.py" --input "$INPUT_DIR" --phagesdb "$PDB_FASTA" --table "$PDB_TABLE" \\', ' --out "$OUT_DIR" --pattern "$PATTERN" --threads "$THREADS" --top "$TOP" --min-contig "$MIN_CONTIG" \\', ' --ncbi "$NCBI_MODE" --ncbi-hosts "$NCBI_HOSTS" --ncbi-fasta "$NCBI_FASTA" $( [[ "$NCBI_REFRESH" == "true" ]] && echo --ncbi-refresh ) \\', ' ${HOST_ARGS[@]+"${HOST_ARGS[@]}"} ${CHECKV_ARGS[@]+"${CHECKV_ARGS[@]}"} ${SIM_ARGS[@]+"${SIM_ARGS[@]}"} \\', ' ${ISO_ARGS[@]+"${ISO_ARGS[@]}"} ${CHECK_ARGS[@]+"${CHECK_ARGS[@]}"}', 'if [[ "$CHECK_ONLY" == "true" ]]; then', ' rm -f "$OUT_DIR/isolation_hosts_pasted.txt"; rmdir "$OUT_DIR" 2>/dev/null || true', ' echo "Run without --check to search."; exit 0', "fi", "", 'REPORT="$OUT_DIR/closest_phages.html"', 'if [[ "$OPEN_REPORT" == "true" && -f "$REPORT" ]] && command -v wslpath >/dev/null 2>&1; then', ' explorer.exe "$(wslpath -w "$REPORT")" >/dev/null 2>&1 || true', 'elif [[ "$OPEN_REPORT" == "true" && -f "$REPORT" && "$(uname)" == "Darwin" ]]; then', ' open "$REPORT" >/dev/null 2>&1 || true', "fi", 'echo "✅ Done."', ")", "", "STATUS=$?", 'if [[ "$STATUS" -ne 0 ]]; then echo; echo "The search stopped with exit code $STATUS. Your terminal stays open."; fi', "" ].join("\n"); } // Downloads folder of the Windows user, from the paths entered (/mnt/c/Users//...) function runCommand(s, extra) { if (isMac()) return `f="$(ls -t ~/Downloads/closest_phages*.sh | head -n 1)"; echo "Running: $f"; caffeinate -i bash "$f"${extra || ""}`; const m = (s.INPUT_DIR || s.PDB_FASTA || "").match(/^(\/mnt\/[a-z]\/Users\/[^/]+)\//); const dl = m ? `${m[1]}/Downloads` : "/mnt/c/Users//Downloads"; return `f="$(ls -t "${dl}"/closest_phages*.sh | head -n 1)"; echo "Running: $f"; bash "$f"${extra || ""}`; } function update() { applyPlatformText(); try { localStorage.setItem("phagePlatform", $("cpPlatform").value); } catch (_) {} for (const [field, preview] of [["inputDir", "inputDirPreview"], ["pdbFasta", "pdbFastaPreview"], ["outDir", "outDirPreview"], ["ncbiFasta", "ncbiFastaPreview"], ["isoHostsFile", "isoHostsFilePreview"]]) { const v = $(field).value.trim(); $(preview).textContent = v ? `${isMac() ? "Path" : "WSL"}: ${toPath(v)}` : (field === "outDir" || field === "ncbiFasta" || field === "isoHostsFile" ? "" : "No path entered"); } const s = readSettings(); const { errors, warnings } = validate(s); const box = $("warning"); box.classList.toggle("cp-error", errors.length > 0); box.style.display = errors.length || warnings.length ? "block" : "none"; box.textContent = [...errors.map((x) => `❌ ${x}`), ...warnings.map((x) => `⚠️ ${x}`)].join("\n"); const ok = errors.length === 0; const script = ok ? makeScript(s) : ""; $("code").textContent = ok ? script : "Enter the genome folder and the PhagesDB FASTA to generate the script."; if (downloadUrl) { URL.revokeObjectURL(downloadUrl); downloadUrl = null; } if (ok) { downloadUrl = URL.createObjectURL(new Blob([script], { type: "text/x-shellscript" })); $("downloadBtn").href = downloadUrl; } else $("downloadBtn").removeAttribute("href"); $("downloadBtn").style.pointerEvents = ok ? "auto" : "none"; $("downloadBtn").style.opacity = ok ? "1" : "0.5"; $("copyBtn").disabled = !ok; $("copyRunBtn").disabled = !ok; $("copyCheckBtn").disabled = !ok; $("copyStatus").textContent = ""; $("copyRunStatus").textContent = ""; { const n = wrap.querySelectorAll(".cp-genus:checked").length; $("genusCount").textContent = `(${n} ticked)`; } $("runCmd").textContent = ok ? runCommand(s) : "Enter the genome folder first."; wrap.value = script; } wrap.querySelectorAll("input, select, textarea").forEach((el) => { el.addEventListener("input", update); el.addEventListener("change", update); }); $("copyRunBtn").addEventListener("click", async () => { try { await navigator.clipboard.writeText($("runCmd").textContent); $("copyRunStatus").textContent = "Copied. Download the script first, then paste the command into " + (isMac() ? "Terminal." : "WSL."); } catch (_) { $("copyRunStatus").textContent = "Copy failed; select the command manually."; } }); $("copyCheckBtn").addEventListener("click", async () => { const s = readSettings(); try { await navigator.clipboard.writeText(runCommand(s, " --check")); $("copyRunStatus").textContent = "Check command copied. Download the script first, then paste it into " + (isMac() ? "Terminal." : "WSL."); } catch (_) { $("copyRunStatus").textContent = "Copy failed; add --check to the end of the run command."; } }); $("copyBtn").addEventListener("click", async () => { try { await navigator.clipboard.writeText($("code").textContent); $("copyStatus").textContent = "Copied."; } catch (_) { $("copyStatus").textContent = "Copy failed; select the script manually."; } }); // ---- v336: remembered form values + Lab settings --------------------------------------------------- { const FORM_KEY = "rudnerForm.closest-phages.v1", LAB_KEY = "rudnerLabSettings.v1"; const LAB_MAP = {"pdbFasta": "phagesdbFasta", "pdbTable": "phagesdbTable", "condaBase": "condaBase", "condaEnv": "medakaEnv", "hostGenome1": "hostGenome1", "hostGenome2": "hostGenome2", "hostGenome3": "hostGenome3", "hostGenome4": "hostGenome4", "hostGenome5": "hostGenome5", "checkvEnv": "checkvEnv", "checkvDb": "checkvDb"}; const GENUS_SEL = ".cp-genus"; const lsGet = k => { try { return JSON.parse(localStorage.getItem(k) || "null"); } catch (_) { return null; } }; const lsSet = (k, v) => { try { localStorage.setItem(k, JSON.stringify(v)); return true; } catch (_) { return false; } }; const fields = () => [...wrap.querySelectorAll("input[id], select[id], textarea[id]")].filter(el => el.type !== "file" && el.type !== "button"); const valOf = el => (el.type === "checkbox" ? el.checked : el.value); const setVal = (el, v) => { if (el.type === "checkbox") el.checked = !!v; else el.value = v; }; const genera = () => (GENUS_SEL ? [...wrap.querySelectorAll(GENUS_SEL + ":checked")].map(x => x.value).sort().join(" ") : ""); const setGenera = list => { const s = new Set(String(list).split(" ")); wrap.querySelectorAll(GENUS_SEL).forEach(x => { x.checked = s.has(x.value); }); }; const defaults = Object.fromEntries(fields().map(el => [el.id, valOf(el)])); const defGenera = genera(); function applyLab(force) { const lab = lsGet(LAB_KEY) || {}, saved = lsGet(FORM_KEY) || {}; let n = 0; for (const [id, key] of Object.entries(LAB_MAP)) { const el = wrap.querySelector("#" + id); if (!el || lab[key] == null || lab[key] === "") continue; if (force || !(id in saved)) { setVal(el, lab[key]); n++; } } return n; } function restore() { const saved = lsGet(FORM_KEY) || {}; for (const el of fields()) if (el.id in saved) setVal(el, saved[el.id]); if (GENUS_SEL && "__genera" in saved) setGenera(saved.__genera); } let timer = null; function save() { clearTimeout(timer); timer = setTimeout(() => { const lab = lsGet(LAB_KEY) || {}, out = {}; for (const el of fields()) { const v = valOf(el), labv = LAB_MAP[el.id] ? lab[LAB_MAP[el.id]] : undefined; if (v !== defaults[el.id] && !(labv != null && labv !== "" && v === labv)) out[el.id] = v; } if (GENUS_SEL && genera() !== defGenera) out.__genera = genera(); lsSet(FORM_KEY, out); }, 300); } const bar = document.createElement("div"); bar.style.cssText = "display:flex;gap:10px;align-items:center;flex-wrap:wrap;font-size:12.5px;color:#5f6b72;margin:0 0 6px"; bar.innerHTML = `Your changes here are remembered in this browser; empty fields come from the Lab settings cell. `; wrap.prepend(bar); const st = bar.querySelector("[data-mem-status]"); bar.querySelector("[data-mem-lab]").addEventListener("click", () => { const n = applyLab(true); st.textContent = n ? `${n} field(s) set from Lab settings.` : "No Lab settings saved yet (add the Lab settings cell to this page)."; save(); update(); }); bar.querySelector("[data-mem-forget]").addEventListener("click", () => { lsSet(FORM_KEY, {}); for (const el of fields()) setVal(el, defaults[el.id]); if (GENUS_SEL) setGenera(defGenera); applyLab(true); st.textContent = "Back to the defaults and Lab settings."; update(); }); applyLab(false); restore(); wrap.addEventListener("input", save); wrap.addEventListener("change", save); window.addEventListener("rudner-lab-settings", () => { applyLab(false); update(); }); } update(); return wrap; }4e:[] 49:{"id":"47b07f65c95574ee","slug":"closest-phages-and-genome-checks","title":"Closest phages and genome checks","thumbnail":"e4a9df7bed4ee0ffca2d57c073f3969b22fbd4a027ab284f2bddb1d6d268d201","access_level":"public","create_time":"$D2026-11-02T18:21:54.109Z","update_time":"$D2026-11-06T20:29:27.294Z","view_time":null,"first_public_time":"$D2026-11-06T19:05:55.481Z","has_unpublished_edits":null,"fork":false,"stars":0,"starred":false,"owner":{"type":"team","id":"d1dfc5f860fb8699","login":"phage","name":"WabushOnce","avatar_url":"https://avatars-next.observableusercontent.com/9303f26b98fd8f7e04b66a6ac3218bfd2ca7e146f15ba24ca32ffac15f802763","tier":"starter_2024","bio":"","followed_by_you":false,"follows_you":false},"version":31,"publish_version":31,"role":null,"authored":false,"authors":[{"type":"individual","id":"3035ce62f81ebcd9","login":"phage","name":"WabushOnce","avatar_url":"https://avatars-next.observableusercontent.com/9303f26b98fd8f7e04b66a6ac3218bfd2ca7e146f15ba24ca32ffac15f802763","tier":null,"bio":"","followed_by_you":false,"follows_you":false}],"members":[{"type":"team","id":"d1dfc5f860fb8699","login":"phage","name":"WabushOnce","avatar_url":"https://avatars-next.observableusercontent.com/9303f26b98fd8f7e04b66a6ac3218bfd2ca7e146f15ba24ca32ffac15f802763","tier":"starter_2024","role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":0,"collection":null,"type":"notebook","body":{"titling":"auto","title":"Closest phages and genome checks","cells":[{"id":1,"value":"# Closest phages and genome checks","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":2,"value":"$4d","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W4e","stdlib":"2"}} 44:null