# Galaxy user-defined tool wrapper for LexicMapStreamer. # # Wraps python/lexicmap_streamer.py from the private repository # https://github.com/nekrut/disassembler (MIT license, (c) Anton Nekrutenko). # Originally vendored verbatim at commit a51eb58d6c1ca941d3fb8d6adf5e8160c3926b91 (repo HEAD # a3999eb60ae78a804d0200b4a1241038077b6e6a, 2026-09-17). v1.0.1 carries one LOCAL correctness # fix on top of that commit, not yet upstreamed to nekrut/disassembler (patch prepared # separately at scratchpad/lexicmap_streamer_fix.diff for the repo owner's own review -- # this is a defect in the user's own private script, not a Foundry- or Galaxy-owned asset, # so it is not filed in this run's foundry-feedback.ledger.yml): # stream_lexicmap_to_msa() no longer assumes LexicMap's tabular output is grouped by accession # (`sgenome`) throughout the whole stream. Real production-scale LexicMap search output is # ranked/ordered by match quality across all matched genomes, not grouped by target genome -- # an accession with more than one HSP/cluster can have its rows scattered non-contiguously # throughout the file. The prior implementation processed one contiguous accession "block" at # a time and raised ValueError on a reappearing accession (confirmed live on real production # data: usegalaxy.org invocation 802eff260023dc62, both lexicmap_streamer_tiling_qc jobs failed # this way against real GENOMIC_PHG/Viral-index results). The fix buffers every row by # accession as it is read and only classifies each accession once the full input has been # consumed, so it is correct regardless of row order -- no upstream pre-sorting required. # Verified: (1) no regression against the original synthetic 3-decoy fixture # (test-data/synthetic_am3/*.tsv, same clean/flagged classification as before); (2) the real # lexicmap_search output that crashed the old script (both gene E and J, ~214k/~215k real # accessions each) now completes in ~7-11s with sane, non-degenerate cohort breakdowns; (3) the # vendored test suite (test_lexicmap_streamer.py, updated: the old # test_ungrouped_input_is_rejected assertion encoded the very assumption being fixed and was # replaced with tests asserting correct merging instead of a raise) passes, 8/8. # # No Tool Shed wrapper exists for this tool (confirmed miss, see galaxy-tool-pin.json); # this UDT is the authoritative wrapper, not a placeholder pending a shed release. # # The script is pure Python 3 standard library (os, sys, gzip, math, json, time, argparse, # urllib.request, pathlib, collections) -- no third-party imports, confirmed by direct read # of both the script and the repo's environment.yml (whose python=3.13 pin and third-party # packages like biopython/requests are for OTHER scripts in the repo, not this one). # # v1.0.2 carries a SECOND local correctness fix, also not yet upstreamed (same rationale as # above -- a defect in the user's own private script). stream_lexicmap_to_msa() now restricts # hit rows to the LexicMap `query` column matching the --ref record, instead of merging every # row in the table onto the reference. LexicMap's `query` tool port is `multiple: true`, so a # whole gene panel connected to it is REDUCED into a single search whose one output table holds # hits for every gene; the previous implementation never read the `query` column and therefore # tiled foreign genes' HSPs onto whichever reference it was given. Confirmed live on # usegalaxy.org history bbd44e69cb8906b54be4b938868d5a9d (invocation 8d8af9e4b30d11f18054bc24110164aa): # all 11 per-gene cohorts had Gene A sequence, truncated to the target gene's length, as their # rank-1 haplotype at 27-33% frequency and ~100% reported pident, inflating Shannon entropy # ~4x (0.38-2.85 -> 4.5-11.9 bits) and suppressing the Gene E am3 `gpE_W7*` signal from 92.59% # to 6.31%. A table with rows but none for --ref now raises instead of emitting empty cohorts. # # v1.0.3 is a MEMORY fix that reverts v1.0.1's mechanism, which the v1.0.2 wiring made both # unnecessary and harmful. v1.0.1 buffered every row by accession and classified only after the # whole table was read -- O(n) in table size. With the workflow now searching ONE gene per job # (nest_gene_panel_for_search) the tables are far larger per gene, and that buffer OOM-killed real # jobs: usegalaxy.org invocation 6a788cc3697550cf, gene Astar, 7.18 GB table, exit 137 with # "Detected 1 oom-kill event(s)"; 3.43 GB passed, so gene A's 9.82 GB table could not have # completed either. v1.0.1's buffering existed only because the pre-Option-A wiring fed one table # containing all 11 genes, where an accession legitimately reappears once per gene. With one gene # per search an accession's HSPs are contiguous again, so streaming is correct: verified on this # run's own output, a 400 MB slice of the real gene C table (444,160 rows, 444,158 distinct # accessions, query == phiX174_C throughout) has ZERO non-contiguous reappearances. v1.0.3 streams # one contiguous accession block at a time in O(1) memory and RAISES on a reappearing accession so # a future violation fails loudly rather than silently double-counting. class: GalaxyUserTool id: lexicmap_streamer version: 1.0.3 name: LexicMap Streamer Multi-HSP Tiling and QC description: >- Streaming multi-HSP coordinate tiling, in-frame codon QC (premature stop-codon screening), 4-tier cohort partitioning (Clean Full-Length / Clean Partial / Flagged Premature Stops / Low-Coverage Fragment), and online haplotype collapsing over a LexicMap tabular search-hit table for one reference coding gene. license: MIT container: python:3.13-slim help: format: markdown content: | Wraps `lexicmap_streamer.py` ("LexicMapStreamer") from [nekrut/disassembler](https://github.com/nekrut/disassembler) (private repository; vendored at commit `a51eb58`, MIT license). For one reference CDS gene and a LexicMap tabular search-hit table, this tool: 0. Selects only the rows whose `query` column equals the reference FASTA's record id (the first token of its header). A LexicMap search submitted with several query sequences returns all of their hits in one table, so this filter is what makes the run gene-specific. Rows for other queries are counted and reported, not merged. If the table has rows but none for this reference, the tool raises rather than emit an empty cohort that would read downstream as a genuine "no hits" result. 1. Groups the selected hit rows by accession and tiles all qualifying HSPs onto the reference coordinate frame, resolving overlaps by the highest-`pident` HSP. 2. Translates the merged sequence in-frame (NCBI code 11) and screens the coding body (excluding the terminal stop) for premature internal stop codons. 3. Partitions each accession into one of four tiers: Clean Full-Length (coverage >= `min_coverage`, 0 stops), Clean Partial (`min_coverage_partial` <= coverage < `min_coverage`, 0 stops), Flagged Premature Stops, or Low-Coverage Fragment. 4. Collapses Clean Full-Length sequences into unique haplotypes as it streams, calling nucleotide/amino-acid mutations against the reference and classifying each haplotype (WT / Synonymous / Missense / Nonsense / Complex). 5. Emits per-gene quantile diagnostics (coverage, identity) and flags a "divergent regime" (median identity < 85% or median coverage < 0.60) recommending de novo graph path-walking (`logan-walker`, out of scope for this workflow) instead. Input requires exactly one reference sequence in the reference FASTA, and a hit table in which all HSP rows for a given accession are contiguous (true of LexicMap output for a single-gene search; the tool raises rather than silently mis-tile a non-contiguous stream). Rows belonging to other queries are filtered out by the `query` column before grouping. Classification streams one accession at a time, holding only that accession's HSPs, so peak memory tracks distinct haplotypes and accession count rather than table size (measured: 132.7 MB vs 527.6 MB on a 400 MB table). `.a2m` output is intentionally suppressed (`--no-a2m`): this workflow declares no A2M port. inputs: - name: reference_fasta type: data format: [fasta] label: Reference CDS FASTA (single gene, in-frame) help: >- The same per-gene reference coding-sequence FASTA used as the LexicMap search query for this gene (one element of the gene query panel). - name: lexicmap_results type: data format: [tabular] label: LexicMap search hits covering this gene help: >- LexicMap's tabular search output (columns must include query, sgenome, pident, qcovHSP, qstart, qend, qseq, sseq). The table may contain hits for several query sequences: only rows whose `query` matches the reference FASTA's record id are used. Rows for a given accession must be contiguous, as LexicMap emits them for a single-gene search. - name: min_coverage type: float value: 0.80 label: Minimum coverage for Clean Full-Length cohort help: Minimum merged reference coverage fraction (0-1) required for the full-length clean tier. - name: min_coverage_partial type: float value: 0.50 label: Minimum coverage for Clean Partial cohort help: Minimum merged reference coverage fraction (0-1) required for the clean partial tier. - name: min_pident type: float value: 60.0 label: Minimum percent identity for candidate HSPs help: HSPs below this percent identity are discarded before tiling. - name: max_internal_stops type: integer label: Maximum premature internal stop codons allowed help: >- Accessions with more premature stop codons than this are excluded from the clean cohorts. No UI default is set here (the confirmed production default is 0, i.e. zero premature stops allowed) because Galaxy's dynamic-tool creation endpoint rejects an integer input whose `value` is exactly `0`; this workflow always supplies the value explicitly via its own `tiling_qc_max_internal_stops` input, so this has no effect on workflow-driven runs. - name: sample_cap type: integer value: 10 label: Max representative sample accessions per haplotype help: Caps how many representative accession IDs are recorded per haplotype in the haplotype TSV. - name: allow_frameshifts type: boolean value: false label: Allow premature-stop/frameshift accessions into the clean cohort help: >- Off by default (confirmed from the source script's argparse default: a store_true flag with no default=True). When off, any accession with more than max_internal_stops premature stops is always routed to the Flagged cohort regardless of coverage. configfiles: - filename: lexicmap_streamer.py content: | #!/usr/bin/env python3 """lexicmap_streamer.py: Streaming multi-HSP coordinate tiling, codon QC, and haplotype collapsing. Streams LexicMap tabular search results directly from usegalaxy.org (over HTTP API), a local TSV file, or standard input (stdin). Core capabilities: 1. Multi-HSP coordinate tiling: Solves the Logan cDBG unitig-split challenge by grouping HSPs by SRA accession (`sgenome`), projecting alignment blocks onto the reference CDS coordinate system, and resolving overlapping positions by highest percent identity. 2. Reading-frame codon QC & Stop-codon screening: Translates in-frame CDS (NCBI Code 11 or Standard), detects premature internal stop codons (e.g. Sanger am3 amber mutants, sequencing errors, or out-of-frame projections), and audits exact codon coordinates. 3. Four-tier cohort partitioning: - Tier 1: Clean Full-Length (coverage >= min_coverage, 0 premature stops) - Tier 2: Clean Partial (min_coverage_partial <= coverage < min_coverage, 0 stops) - Tier 3: Flagged Premature Stops (stops > max_stops) - Tier 4: Low-Coverage Fragments (coverage < min_coverage_partial) 4. Online haplotype collapser: Collapses clean full-length sequences in O(1) streaming memory, calls nucleotide and amino acid mutations, classifies mutation impact (WT, Synonymous, Missense, Nonsense), and outputs clean MSAs, mutation frequency tables, and complete cohort ledgers. 5. Regime diagnostics: Calculates coverage and identity quantiles (p10, p25, median, p75, p90), detects divergent targets (median pid < 85% or cov < 0.60), and alerts when de novo graph path-walking with `logan-walker` is required. Outputs generated: - .clean.msa.fasta (collapsed clean haplotypes with rank and count) - .clean.haplotypes.tsv (detailed mutation audit, frequencies, sample accessions) - .clean.accessions.fasta (uncollapsed full-length clean accessions) - .clean_expanded.accessions.fasta (clean cohort including >=50% partials) - .flagged.accessions.fasta (flagged sequences: stops, frameshifts, or low cov) - .flagged.tsv (audit table of all flagged accessions and reasons) - .cohort_ledger.tsv (complete audit ledger of all evaluated accessions) - .summary.json (quantiles, cohort breakdown, multi-HSP rescue counts) - .clean.a2m (A2M formatted MSA for protein language models / EVcouplings) Examples: # Stream directly from a Galaxy dataset ID over HTTP: python/lexicmap_streamer.py \\ -q queries/markers/zika_ns5.fasta \\ -d f9cad7b01a47213500416a19986dca5c \\ -o out/zika_lexicmap # Stream from a local LexicMap tabular file: python/lexicmap_streamer.py \\ -q queries/query_hiv1_pol.fasta \\ -f /data/hiv1_pol_lexicmap.tsv \\ -o out/hiv1_pol # Pipe directly from a compressed file or stream via STDIN: zcat /data/lexicmap_output.tsv.gz | python/lexicmap_streamer.py \\ -q queries/master/pfcrt.fasta \\ -f - \\ -o out/pfcrt """ import os import sys import gzip import math import json import time import argparse import urllib.request from pathlib import Path from collections import Counter, defaultdict # Translation table: NCBI standard / bacterial code 11 CODON_TABLE_11 = { 'TTT': 'F', 'TTC': 'F', 'TTA': 'L', 'TTG': 'L', 'TCT': 'S', 'TCC': 'S', 'TCA': 'S', 'TCG': 'S', 'TAT': 'Y', 'TAC': 'Y', 'TAA': '*', 'TAG': '*', 'TGT': 'C', 'TGC': 'C', 'TGA': '*', 'TGG': 'W', 'CTT': 'L', 'CTC': 'L', 'CTA': 'L', 'CTG': 'L', 'CCT': 'P', 'CCC': 'P', 'CCA': 'P', 'CCG': 'P', 'CAT': 'H', 'CAC': 'H', 'CAA': 'Q', 'CAG': 'Q', 'CGT': 'R', 'CGC': 'R', 'CGA': 'R', 'CGG': 'R', 'ATT': 'I', 'ATC': 'I', 'ATA': 'I', 'ATG': 'M', 'ACT': 'T', 'ACC': 'T', 'ACA': 'T', 'ACG': 'T', 'AAT': 'N', 'AAC': 'N', 'AAA': 'K', 'AAG': 'K', 'AGT': 'S', 'AGC': 'S', 'AGA': 'R', 'AGG': 'R', 'GTT': 'V', 'GTC': 'V', 'GTA': 'V', 'GTG': 'V', 'GCT': 'A', 'GCC': 'A', 'GCA': 'A', 'GCG': 'A', 'GAT': 'D', 'GAC': 'D', 'GAA': 'E', 'GAG': 'E', 'GGT': 'G', 'GGC': 'G', 'GGA': 'G', 'GGG': 'G' } def get_galaxy_api_key(api_key: str = None, key_file: str = None) -> str: """Resolves Galaxy API key from CLI, environment, or default credentials file.""" if api_key: return api_key.strip() if os.environ.get("GALAXY_API_KEY"): return os.environ["GALAXY_API_KEY"].strip() if key_file and Path(key_file).exists(): return Path(key_file).read_text().strip() default_tmp = Path("/tmp/gxy.txt") if default_tmp.exists(): return default_tmp.read_text().strip() raise RuntimeError( "Galaxy API key not found. Provide --api-key, set GALAXY_API_KEY, or write /tmp/gxy.txt." ) def translate_nt(seq: str, codon_table: dict = CODON_TABLE_11) -> str: """Translates nucleotide sequence to amino acids.""" seq = seq.upper() aa_list = [] for i in range(0, len(seq) - 2, 3): codon = seq[i:i + 3] if '-' in codon or len(codon) < 3: aa_list.append('-') else: aa_list.append(codon_table.get(codon, 'X')) return "".join(aa_list) def find_internal_stops(seq: str, ref_len: int, codon_table: dict = CODON_TABLE_11) -> list: """ Finds 1-based codon positions of premature stop codons in coding sequence body. Excludes the terminal natural stop codon (last 3 nt) of the CDS. """ stop_positions = [] for i in range(0, ref_len - 3, 3): codon = seq[i:i + 3] if '-' not in codon and len(codon) == 3: if codon_table.get(codon, 'X') == '*': stop_positions.append(i // 3 + 1) return stop_positions def compute_quantiles(values: list) -> dict: """Calculates summary percentiles for a metric.""" if not values: return {"p10": 0.0, "p25": 0.0, "p50_median": 0.0, "p75": 0.0, "p90": 0.0, "mean": 0.0} s = sorted(values) n = len(s) def q(p): idx = int(round(p * (n - 1))) return round(s[max(0, min(n - 1, idx))], 4) return { "p10": q(0.10), "p25": q(0.25), "p50_median": q(0.50), "p75": q(0.75), "p90": q(0.90), "mean": round(sum(s) / n, 4) } def load_reference(fasta_path: str): """Loads reference sequence, length, and header from FASTA file.""" path = Path(fasta_path) if not path.exists(): raise FileNotFoundError(f"Reference FASTA file not found: {path}") header = "" seq_parts = [] with open(path, "r") as f: for line in f: line = line.strip() if not line: continue if line.startswith(">"): if not header: header = line[1:] else: seq_parts.append(line) ref_seq = "".join(seq_parts).upper() ref_name = header.split()[0] if header else path.stem return ref_name, header, ref_seq def get_galaxy_stream(dataset_id: str, base_url: str = "https://usegalaxy.org", api_key: str = None, key_file: str = None): """Opens a buffered HTTP streaming connection to Galaxy dataset display endpoint.""" key = get_galaxy_api_key(api_key, key_file) url = f"{base_url.rstrip('/')}/api/datasets/{dataset_id}/display" req = urllib.request.Request(url, headers={"x-api-key": key}) response = urllib.request.urlopen(req, timeout=120) return response def resolve_dataset_state(dataset_id: str, base_url: str = "https://usegalaxy.org", api_key: str = None, key_file: str = None): """Checks the live state of a Galaxy dataset.""" try: key = get_galaxy_api_key(api_key, key_file) url = f"{base_url.rstrip('/')}/api/datasets/{dataset_id}" req = urllib.request.Request(url, headers={"x-api-key": key}) with urllib.request.urlopen(req, timeout=30) as res: data = json.loads(res.read()) return data.get("state") except Exception: return None def merge_accession_hsps(hsps: list, ref_len: int, min_pident: float, min_coverage: float): """ Merges all qualifying HSPs of an accession onto the reference coordinate system. Overlapping positions between HSPs are resolved by selecting the base from the HSP with higher pident. """ grid = [None] * ref_len single_ge_min_cov = False single_full = False for pident, qcov, qstart, qend, qseq, sseq in hsps: if pident < min_pident: continue if qcov >= (min_coverage * 100.0): single_ge_min_cov = True if qstart == 1 and qend == ref_len and pident >= min_pident: single_full = True q_pos = qstart - 1 for qc, sc in zip(qseq, sseq): if qc != '-': if 0 <= q_pos < ref_len: ex = grid[q_pos] if ex is None or pident > ex[1]: grid[q_pos] = (sc.upper(), pident) q_pos += 1 # Mean identity over the positions that actually carry a subject base. The numerator must # use the same non-gap positions as `covered_bases`; summing pident over gap positions too # (subject deletions, where x[0] == '-') inflates the mean and can push it above 100%. covered = [x for x in grid if x is not None and x[0] != '-'] covered_bases = len(covered) mean_pid = (sum(x[1] for x in covered) / covered_bases) if covered_bases > 0 else 0.0 merged_seq = "".join(grid[i][0] if grid[i] is not None else '-' for i in range(ref_len)) return merged_seq, covered_bases, mean_pid, single_ge_min_cov, single_full def stream_lexicmap_to_msa( ref_path: str, dataset_id: str = None, file_path: str = None, base_url: str = "https://usegalaxy.org", api_key: str = None, api_key_file: str = "/tmp/gxy.txt", min_coverage: float = 0.80, min_coverage_partial: float = 0.50, min_pident: float = 60.0, max_internal_stops: int = 0, allow_frameshifts: bool = False, full_length_only: bool = False, max_records: int = None, sample_cap: int = 10, out_prefix: str = None, write_a2m: bool = True ) -> dict: """ Main streaming pipeline: - Streams rows, groups by accession, merges multi-HSP alignments - Partitions cohort into: 1. Clean Full-Length (cov >= min_coverage, 0 internal stops) 2. Clean Partial (min_coverage_partial <= cov < min_coverage, 0 stops) 3. Flagged (internal stops, frameshifts, or low coverage) - Writes clean MSAs, flagged audit FASTA/TSV, and complete cohort ledger. """ ref_name, ref_header, ref_seq = load_reference(ref_path) ref_len = len(ref_seq) ref_is_cds = (ref_len % 3 == 0) ref_aa = translate_nt(ref_seq) if ref_is_cds else "" if dataset_id: ds_state = resolve_dataset_state(dataset_id, base_url, api_key, api_key_file) if ds_state and ds_state != "ok": print(f"[!] Warning: Galaxy dataset {dataset_id} is in state '{ds_state}' (not 'ok').") if ds_state == "queued": print(" Job is still queued on Galaxy. Skipping for now.") return {"status": "queued", "dataset_id": dataset_id} print("=" * 80) print(" LexicMapStreamer: Multi-HSP Coordinate Tiler & Cohort Partitioner") print("=" * 80) print(f" Reference Sequence : {ref_name} ({ref_len} bp, {len(ref_aa)} codons)") print(f" Reference Header : {ref_header}") print(f" Target QC Settings : full_length_cov={min_coverage*100:.1f}%, partial_cov={min_coverage_partial*100:.1f}%, " f"min_pident={min_pident:.1f}%, max_stops={max_internal_stops if ref_is_cds else 'N/A'}") if dataset_id: print(f" Input Source : Streaming Galaxy dataset: {dataset_id} ({base_url})") stream = get_galaxy_stream(dataset_id, base_url, api_key, api_key_file) is_http = True elif file_path == "-": print(" Input Source : Reading from STDIN") stream = sys.stdin is_http = False elif file_path: print(f" Input Source : Local file: {file_path}") if file_path.endswith(".gz"): stream = gzip.open(file_path, "rt", encoding="utf-8") else: stream = open(file_path, "r", encoding="utf-8") is_http = False else: raise ValueError("Must specify either --dataset or --file.") if not out_prefix: out_dir = Path("msa") out_dir.mkdir(parents=True, exist_ok=True) out_prefix = f"msa/{ref_name}" else: out_path = Path(out_prefix) if out_path.is_dir() or str(out_prefix).endswith("/"): out_path.mkdir(parents=True, exist_ok=True) out_prefix = str(out_path / ref_name) else: out_path.parent.mkdir(parents=True, exist_ok=True) # Read header line (skipping comment lines if any) while True: header_line = stream.readline() if not header_line: raise ValueError("Input stream was empty or did not contain a header line.") if is_http and isinstance(header_line, bytes): header_line = header_line.decode("utf-8") header_line = header_line.strip() if header_line and not (header_line.startswith("#") and "\t" not in header_line): break # Strip leading '#' if present on first column if header_line.startswith("#"): header_line = header_line[1:].strip() header_cols = header_line.split("\t") col_idx = {c.strip(): i for i, c in enumerate(header_cols)} # Schema resolution def find_col(candidates): for c in candidates: if c in col_idx: return col_idx[c] # Case-insensitive fallback lower_map = {k.lower(): v for k, v in col_idx.items()} for c in candidates: if c.lower() in lower_map: return lower_map[c.lower()] return None idx_query = find_col(["query", "qseqid", "qid", "query_id"]) idx_sgenome = find_col(["sgenome", "subject", "sseqid", "acc", "accession"]) idx_pident = find_col(["pident", "pid", "perc_identity", "identity"]) idx_qcov = find_col(["qcovHSP", "qcov", "qcovhsp", "q_cov"]) idx_qstart = find_col(["qstart", "q_start", "query_start"]) idx_qend = find_col(["qend", "q_end", "query_end"]) idx_qseq = find_col(["qseq", "query_seq", "q_seq"]) idx_sseq = find_col(["sseq", "subject_seq", "s_seq"]) missing = [] if idx_query is None: missing.append("query") if idx_sgenome is None: missing.append("sgenome/subject") if idx_pident is None: missing.append("pident") if idx_qcov is None: missing.append("qcovHSP") if idx_qstart is None: missing.append("qstart") if idx_qend is None: missing.append("qend") if idx_qseq is None: missing.append("qseq") if idx_sseq is None: missing.append("sseq") if missing: raise KeyError(f"Missing required columns {missing} in LexicMap header: {header_cols}") # Cohort file handles f_clean_fa = open(f"{out_prefix}.clean.accessions.fasta", "w", encoding="utf-8") f_clean_exp_fa = open(f"{out_prefix}.clean_expanded.accessions.fasta", "w", encoding="utf-8") f_flagged_fa = open(f"{out_prefix}.flagged.accessions.fasta", "w", encoding="utf-8") f_flagged_tsv = open(f"{out_prefix}.flagged.tsv", "w", encoding="utf-8") f_cohort_tsv = open(f"{out_prefix}.cohort_ledger.tsv", "w", encoding="utf-8") f_flagged_tsv.write("accession\tcoverage\tmean_pident\tflag\tn_stops\tstop_codons\thsps_merged\n") f_cohort_tsv.write("accession\tcoverage\tmean_pident\tstatus\ttier\tn_stops\tstop_codons\thsps_merged\n") # Tracking & State clean_full_haps = Counter() clean_full_pident_sum = defaultdict(float) clean_full_samples = defaultdict(list) all_covs = [] all_pids = [] cov_ge_50_count = 0 cov_ge_50_stops = 0 total_rows = 0 rows_skipped_other_query = 0 queries_seen = {} total_accessions = 0 count_clean_full = 0 count_clean_partial = 0 count_flagged_stops = 0 count_flagged_low_cov = 0 total_rescued_multihsp = 0 total_single_hsp_full = 0 start_time = time.time() last_report_time = start_time print("\n---> Streaming rows, grouping by accession, and partitioning into cohorts...") # MEMORY FIX (v1.0.3): v1.0.1/v1.0.2 buffered EVERY row by accession and only classified # once the whole input had been consumed. That was O(n) in the table size and OOM-killed # real jobs: usegalaxy.org invocation 6a788cc3697550cf, gene Astar, a 7.18 GB hit table, # exit 137 ("Detected 1 oom-kill event(s)"). Tables of 3.43 GB passed and 7.18 GB died, so # gene A at 9.82 GB could not have completed either. # # The buffering existed only because the pre-Option-A wiring fed ONE table containing all # 11 genes' hits, in which an accession legitimately reappears once per gene. Now that # nest_gene_panel_for_search gives each search a single gene, an accession's HSPs are # contiguous again and streaming is correct. Verified on this run's own output before # making the change: a 400 MB slice of the real gene C table (444,160 rows, 444,158 # distinct accessions, query == phiX174_C throughout) contains ZERO non-contiguous # accession reappearances. # # Process one contiguous accession block at a time: memory for the row stream itself is # O(1) (only the current accession's HSPs are held). Note the job is NOT O(1) overall -- # the haplotype counter and the coverage/identity quantile lists still grow with distinct # haplotypes and accession count -- but that is far smaller than holding every row, and it # is the same profile as the reference implementation. Measured on a real 400 MB gene C # table: 527.6 MB peak RSS buffering vs 132.7 MB streaming, byte-identical results. Keep a guard that # RAISES on a reappearing accession: if the precondition is ever violated again the job # fails loudly instead of silently double-counting the accession and undercounting its # multi-HSP coverage. curr_acc = None curr_hsps = [] seen_accessions = set() def process_accession_block(acc, hsps): nonlocal total_accessions, count_clean_full, count_clean_partial nonlocal count_flagged_stops, count_flagged_low_cov nonlocal total_rescued_multihsp, total_single_hsp_full nonlocal cov_ge_50_count, cov_ge_50_stops if not acc or not hsps: return if acc in seen_accessions: raise ValueError( f"Accession '{acc}' reappears after its block was already processed, so this " f"hit table is not grouped by accession. All HSPs of an accession must be " f"contiguous for streaming classification. This should not happen for a " f"single-gene search; if the table covers several query genes, that is the " f"upstream defect to fix. As a stopgap, sort the table by the sgenome column " f"first (`sort -k -s`) and re-run." ) seen_accessions.add(acc) total_accessions += 1 merged_seq, covered_bases, mean_pid, single_ge_cov, single_full = merge_accession_hsps( hsps, ref_len, min_pident, min_coverage ) cov = covered_bases / ref_len all_covs.append(cov) if cov > 0: all_pids.append(mean_pid) # Check premature stop codons in coding body stop_codons = find_internal_stops(merged_seq, ref_len) if ref_is_cds else [] n_stops = len(stop_codons) stop_str = ",".join(str(p) for p in stop_codons) if stop_codons else "None" if cov >= 0.50: cov_ge_50_count += 1 if n_stops > 0: cov_ge_50_stops += 1 # Stratify accession. With --full-length-only, the clean full-length tier requires # complete coverage of the reference (every position tiled), not merely min_coverage. full_cov = 1.0 if full_length_only else min_coverage if cov >= full_cov and (n_stops <= max_internal_stops or allow_frameshifts): status = "CLEAN_FULL_LENGTH" tier = 1 count_clean_full += 1 if single_full: total_single_hsp_full += 1 if not single_ge_cov: total_rescued_multihsp += 1 clean_full_haps[merged_seq] += 1 clean_full_pident_sum[merged_seq] += mean_pid samples = clean_full_samples[merged_seq] if len(samples) < sample_cap and acc not in samples: samples.append(acc) f_clean_fa.write(f">{acc} cov={cov:.4f} pid={mean_pid:.2f} status=CLEAN_FULL_LENGTH\n{merged_seq}\n") f_clean_exp_fa.write(f">{acc} cov={cov:.4f} pid={mean_pid:.2f} status=CLEAN_FULL_LENGTH\n{merged_seq}\n") elif cov >= min_coverage_partial and (n_stops <= max_internal_stops or allow_frameshifts): status = "CLEAN_PARTIAL" tier = 2 count_clean_partial += 1 f_clean_exp_fa.write(f">{acc} cov={cov:.4f} pid={mean_pid:.2f} status=CLEAN_PARTIAL\n{merged_seq}\n") elif n_stops > max_internal_stops and not allow_frameshifts: status = "FLAGGED_PREMATURE_STOPS" tier = 3 count_flagged_stops += 1 f_flagged_fa.write(f">{acc} cov={cov:.4f} pid={mean_pid:.2f} flag=PREMATURE_STOPS stops={n_stops} codons={stop_str}\n{merged_seq}\n") f_flagged_tsv.write(f"{acc}\t{cov:.4f}\t{mean_pid:.2f}\tPREMATURE_STOPS\t{n_stops}\t{stop_str}\t{len(hsps)}\n") else: # cov < min_coverage_partial and no stops status = "LOW_COVERAGE_FRAGMENT" tier = 4 count_flagged_low_cov += 1 f_flagged_fa.write(f">{acc} cov={cov:.4f} pid={mean_pid:.2f} flag=LOW_COVERAGE\n{merged_seq}\n") f_flagged_tsv.write(f"{acc}\t{cov:.4f}\t{mean_pid:.2f}\tLOW_COVERAGE\t0\tNone\t{len(hsps)}\n") f_cohort_tsv.write(f"{acc}\t{cov:.4f}\t{mean_pid:.2f}\t{status}\t{tier}\t{n_stops}\t{stop_str}\t{len(hsps)}\n") try: while True: line = stream.readline() if not line: break if is_http and isinstance(line, bytes): line = line.decode("utf-8") line = line.strip() if not line: continue parts = line.split("\t") if len(parts) <= max(idx_query, idx_sgenome, idx_pident, idx_qcov, idx_qstart, idx_qend, idx_qseq, idx_sseq): continue # CORRECTNESS FIX (v1.0.2): a single LexicMap search may carry hits for MANY # query sequences -- the `query` tool port is `multiple: true`, so a whole # gene panel submitted at once lands in ONE result table. Rows must therefore # be restricted to the query matching THIS run's --ref record before being # merged, or HSPs from unrelated genes are tiled onto this reference's # coordinates and silently corrupt every downstream cohort, haplotype and # entropy figure. (Observed live: usegalaxy.org history # bbd44e69cb8906b54be4b938868d5a9d, where all 11 per-gene cohorts were # dominated by Gene A sequence truncated to the target gene's length, at # ~100% reported pident.) Prior versions ignored the `query` column entirely. row_query = parts[idx_query] queries_seen[row_query] = queries_seen.get(row_query, 0) + 1 if row_query != ref_name: rows_skipped_other_query += 1 continue total_rows += 1 acc = parts[idx_sgenome] pident = float(parts[idx_pident]) qcov = float(parts[idx_qcov]) qstart = int(parts[idx_qstart]) qend = int(parts[idx_qend]) qseq = parts[idx_qseq] sseq = parts[idx_sseq] hsp = (pident, qcov, qstart, qend, qseq, sseq) if acc != curr_acc: if curr_acc is not None: process_accession_block(curr_acc, curr_hsps) curr_acc = acc curr_hsps = [hsp] else: curr_hsps.append(hsp) now = time.time() if now - last_report_time >= 5.0 or (max_records and total_rows % 25000 == 0): elapsed = now - start_time rate = total_rows / elapsed if elapsed > 0 else 0 print(f" [Stream] Rows for {ref_name}: {total_rows:,} | Skipped (other query): " f"{rows_skipped_other_query:,} | Accs: {total_accessions:,} | " f"Clean Full: {count_clean_full:,} | Rescued: {total_rescued_multihsp:,} | " f"Flagged Stops: {count_flagged_stops:,} | Speed: {rate:,.0f} rows/s") last_report_time = now if max_records and total_rows >= max_records: print(f"\n[*] Reached row limit (--max-records {max_records:,}).") break # A result table that contains rows, but none for this reference, is never a # legitimate empty cohort -- it means --ref does not correspond to any query in # this search (renamed FASTA header, wrong hit table wired to the step, or a # panel/reference mismatch). Fail loudly rather than emit empty cohorts that look # like a real "no hits" result downstream. if total_rows == 0 and rows_skipped_other_query > 0: observed = ", ".join( f"{q} ({n:,} rows)" for q, n in sorted(queries_seen.items(), key=lambda kv: -kv[1])[:10] ) raise ValueError( f"No rows in the hit table have query == '{ref_name}' (the first token of " f"the --ref FASTA header), but {rows_skipped_other_query:,} rows for other " f"queries were present. The reference and the hit table do not correspond. " f"Queries actually present: {observed}." ) if rows_skipped_other_query > 0: print(f"\n---> Query filter: kept {total_rows:,} rows for '{ref_name}'; skipped " f"{rows_skipped_other_query:,} rows belonging to {len(queries_seen) - 1} " f"other quer{'y' if len(queries_seen) == 2 else 'ies'} in the same hit table.") # Flush the final accession block. if curr_acc is not None: process_accession_block(curr_acc, curr_hsps) finally: if hasattr(stream, "close"): stream.close() f_clean_fa.close() f_clean_exp_fa.close() f_flagged_fa.close() f_flagged_tsv.close() f_cohort_tsv.close() total_time = time.time() - start_time # Distributions & Diagnostics cov_stats = compute_quantiles(all_covs) pid_stats = compute_quantiles(all_pids) med_cov = cov_stats["p50_median"] med_pid = pid_stats["p50_median"] stop_rate_ge_50 = (cov_ge_50_stops / cov_ge_50_count * 100.0) if cov_ge_50_count > 0 else 0.0 is_divergent = (med_pid < 85.0 or med_cov < 0.60) regime = "divergent" if is_divergent else "conserved" # A valid stream can carry a header and no data rows; guard every by-total percentage so an # empty cohort reports zeros rather than raising ZeroDivisionError. def pct_of_total(n): return (n / total_accessions * 100.0) if total_accessions > 0 else 0.0 print(f"\n[✓] Stream processing complete in {total_time:.2f} seconds.") print(f" Total rows streamed : {total_rows:,}") print(f" Rows skipped (other query) : {rows_skipped_other_query:,}") print(f" Total unique SRA accessions : {total_accessions:,}") print(f" 1. Clean Full-Length Cohort : {count_clean_full:,} accessions ({pct_of_total(count_clean_full):.2f}% of total)") print(f" - Multi-HSP rescued : {total_rescued_multihsp:,} accessions") print(f" - Single-HSP full-length : {total_single_hsp_full:,} accessions") print(f" 2. Clean Partials (50%-80% cov): {count_clean_partial:,} accessions") print(f" --> Combined Clean Cohort : {count_clean_full + count_clean_partial:,} accessions ({pct_of_total(count_clean_full + count_clean_partial):.2f}%)") print(f" 3. Flagged for Premature Stops: {count_flagged_stops:,} accessions") print(f" 4. Low-Coverage Fragments (<50%): {count_flagged_low_cov:,} accessions") print(f" Coverage Distribution (All) : p10={cov_stats['p10']:.2f}, median={med_cov:.2f}, p90={cov_stats['p90']:.2f}") print(f" Identity Distribution (All) : p10={pid_stats['p10']:.1f}%, median={med_pid:.1f}%, p90={pid_stats['p90']:.1f}%") warning_msg = None if is_divergent: warning_msg = ( f"DIVERGENT REGIME DETECTED: Median identity is {med_pid:.1f}% (< 85%) and/or " f"median coverage is {med_cov:.2f} (< 0.60) with {stop_rate_ge_50:.1f}% premature stops in >=50% covered accessions. " f"Reference projection is fragmentary. Clean full-length sequences ({count_clean_full:,}) and clean partials " f"({count_clean_partial:,}) have been safely partitioned from flagged stop-codon projections ({count_flagged_stops:,}). " f"For unbiased whole-cohort haplotype assembly on this target, use disassembler (logan-walker)." ) print("\n" + "!" * 80) print(" [!] DIVERGENT TARGET REGIME DETECTED") print("!" * 80) print(warning_msg) print("!" * 80 + "\n") # Sort clean full-length haplotypes sorted_haplotypes = clean_full_haps.most_common() # Compute Shannon entropy shannon_entropy = 0.0 for seq, count in sorted_haplotypes: p = count / count_clean_full if count_clean_full > 0 else 0 if p > 0: shannon_entropy -= p * math.log2(p) msa_clean_fasta = Path(f"{out_prefix}.clean.msa.fasta") msa_compat_fasta = Path(f"{out_prefix}.msa.fasta") counts_tsv_path = Path(f"{out_prefix}.clean.haplotypes.tsv") counts_compat_tsv = Path(f"{out_prefix}.haplotypes.tsv") summary_json_path = Path(f"{out_prefix}.summary.json") print(f"\n---> Emitting clean alignments and cohort tables to '{out_prefix}'...") # Write clean MSA FASTA for path in [msa_clean_fasta, msa_compat_fasta]: with open(path, "w", encoding="utf-8") as f_fa: f_fa.write(f">REFERENCE_{ref_name} len={ref_len} count={count_clean_full} type=reference\n{ref_seq}\n") for rank, (seq, count) in enumerate(sorted_haplotypes, 1): freq = count / count_clean_full if count_clean_full > 0 else 0 is_wt = (seq == ref_seq) tag = "WT" if is_wt else f"variant_rank_{rank}" f_fa.write(f">haplotype_{rank:05d} count={count} freq={freq:.5f} tag={tag}\n{seq}\n") # Write A2M if write_a2m: with open(f"{out_prefix}.clean.a2m", "w", encoding="utf-8") as f_a2m, open(f"{out_prefix}.a2m", "w", encoding="utf-8") as f_a2m_compat: for f in [f_a2m, f_a2m_compat]: f.write(f">{ref_name}/1-{ref_len}\n{ref_seq}\n") for rank, (seq, count) in enumerate(sorted_haplotypes, 1): f.write(f">hap_{rank}_count_{count}/1-{ref_len}\n{seq}\n") # Write Haplotype TSV top_preview = [] for tsv_p in [counts_tsv_path, counts_compat_tsv]: with open(tsv_p, "w", encoding="utf-8") as f_tsv: headers = [ "haplotype_id", "rank", "count", "frequency", "mean_pident", "coverage_pct", "num_nt_muts", "nt_mutations", "num_aa_muts", "aa_mutations", "mutation_class", "representative_samples", "sequence" ] f_tsv.write("\t".join(headers) + "\n") for rank, (seq, count) in enumerate(sorted_haplotypes, 1): freq = count / count_clean_full if count_clean_full > 0 else 0 mean_pid = clean_full_pident_sum[seq] / count if count > 0 else 0 cov_pct = sum(1 for c in seq if c != '-') / ref_len * 100.0 nt_muts = [] for i, (r_char, h_char) in enumerate(zip(ref_seq, seq)): if h_char != '-' and h_char != r_char: nt_muts.append(f"{i+1}{r_char}>{h_char}") aa_muts = [] mut_class = "WT" if seq == ref_seq: mut_class = "WT" elif ref_is_cds: h_aa = translate_nt(seq) for r_idx, (r_a, h_a) in enumerate(zip(ref_aa, h_aa)): if h_a != '-' and h_a != r_a: aa_muts.append(f"{r_a}{r_idx+1}{h_a}") if not aa_muts and not nt_muts: mut_class = "WT" elif not aa_muts and nt_muts: mut_class = "Synonymous" elif any(a.endswith("*") for a in aa_muts): mut_class = "Nonsense" elif aa_muts: mut_class = "Missense" else: mut_class = "Complex" else: mut_class = "NonCoding_WT" if not nt_muts else "NonCoding_Variant" nt_mut_str = ",".join(nt_muts) if nt_muts else "WT" aa_mut_str = ",".join(aa_muts) if aa_muts else ("WT" if mut_class in ("WT", "Synonymous") else "NA") samples_str = ",".join(clean_full_samples[seq]) row_data = [ f"haplotype_{rank:05d}", str(rank), str(count), f"{freq:.6f}", f"{mean_pid:.2f}", f"{cov_pct:.2f}", str(len(nt_muts)), nt_mut_str, str(len(aa_muts)), aa_mut_str, mut_class, samples_str, seq ] f_tsv.write("\t".join(row_data) + "\n") if rank <= 10 and tsv_p == counts_tsv_path: top_preview.append({ "rank": rank, "count": count, "freq": freq, "mut_class": mut_class, "nt_muts": nt_mut_str, "aa_muts": aa_mut_str, "samples": samples_str }) # Summary JSON summary_data = { "status": "ok", "reference_name": ref_name, "reference_length_bp": ref_len, "is_cds": ref_is_cds, "regime": regime, "regime_warning": warning_msg, "total_streamed_rows": total_rows, "rows_skipped_other_query": rows_skipped_other_query, "total_accessions_evaluated": total_accessions, "cohort_breakdown": { "clean_full_length": count_clean_full, "clean_partial": count_clean_partial, "combined_clean_cohort": count_clean_full + count_clean_partial, "flagged_premature_stops": count_flagged_stops, "flagged_low_coverage": count_flagged_low_cov }, "recovery_metrics": { "multi_hsp_rescued_accessions": total_rescued_multihsp, "single_hsp_full_length_accessions": total_single_hsp_full, "premature_stop_rate_pct_at_50cov": round(stop_rate_ge_50, 2) }, "coverage_quantiles": cov_stats, "pident_quantiles": pid_stats, "clean_unique_haplotypes": len(sorted_haplotypes), "shannon_entropy_bits": round(shannon_entropy, 4), "wildtype_count": clean_full_haps.get(ref_seq, 0), "wildtype_fraction": round(clean_full_haps.get(ref_seq, 0) / count_clean_full, 6) if count_clean_full > 0 else 0, "top_10_haplotypes": top_preview, "runtime_seconds": round(total_time, 2) } summary_json_path.write_text(json.dumps(summary_data, indent=2)) print(f" [✓] Clean MSA FASTA : {msa_clean_fasta} ({len(sorted_haplotypes) + 1} sequences)") print(f" [✓] Clean Full Accessions : {out_prefix}.clean.accessions.fasta ({count_clean_full:,} accessions)") print(f" [✓] Clean Expanded (>=50%): {out_prefix}.clean_expanded.accessions.fasta ({count_clean_full + count_clean_partial:,} accessions)") print(f" [✓] Flagged Accessions : {out_prefix}.flagged.accessions.fasta ({count_flagged_stops + count_flagged_low_cov:,} accessions)") print(f" [✓] Flagged Audit Table : {out_prefix}.flagged.tsv") print(f" [✓] Complete Cohort Ledger: {out_prefix}.cohort_ledger.tsv ({total_accessions:,} rows)") print(f" [✓] Summary JSON : {summary_json_path}") # Display Top Table print("\n" + "=" * 80) print(f" Top Clean Haplotypes ({ref_name})") print("=" * 80) print(f"{'Rank':<5} {'Count':<9} {'Freq':<8} {'Class':<12} {'AA Mutations':<18} {'NT Mutations':<25}") print("-" * 80) for row in top_preview: aa_disp = row['aa_muts'][:16] + ".." if len(row['aa_muts']) > 16 else row['aa_muts'] nt_disp = row['nt_muts'][:23] + ".." if len(row['nt_muts']) > 23 else row['nt_muts'] print(f"{row['rank']:<5} {row['count']:<9} {row['freq']:<8.4f} {row['mut_class']:<12} {aa_disp:<18} {nt_disp:<25}") print("=" * 80) return summary_data def parse_args(): parser = argparse.ArgumentParser( description="LexicMapStreamer: Stream LexicMap search results, tile multi-HSP alignments, stratify cohorts, and collapse haplotypes." ) parser.add_argument("-q", "-r", "--query", "--ref", dest="ref", type=str, required=True, help="Path to reference FASTA file (e.g. queries/query_hiv1_pol.fasta).") parser.add_argument("-d", "--dataset", type=str, default=None, help="Galaxy dataset ID to stream from over HTTP.") parser.add_argument("-f", "-i", "--file", "--input", dest="file", type=str, default=None, help="Path to local LexicMap tabular file (or '-' for stdin).") parser.add_argument("--base-url", type=str, default="https://usegalaxy.org", help="Galaxy server base URL (default: https://usegalaxy.org).") parser.add_argument("--api-key", type=str, default=None, help="Galaxy API key string.") parser.add_argument("--api-key-file", type=str, default="/tmp/gxy.txt", help="Path to file containing Galaxy API key (default: /tmp/gxy.txt).") parser.add_argument("--min-coverage", type=float, default=0.80, help="Minimum merged query coverage fraction for full-length clean cohort (default: 0.80 = 80%%).") parser.add_argument("--min-coverage-partial", type=float, default=0.50, help="Minimum coverage fraction for clean partial expanded cohort (default: 0.50 = 50%%).") parser.add_argument("--min-pident", type=float, default=60.0, help="Minimum percent identity for candidate HSPs (default: 60.0%%).") parser.add_argument("--max-internal-stops", type=int, default=0, help="Maximum premature stop codons allowed in CDS body (default: 0).") parser.add_argument("--allow-frameshifts", action="store_true", help="Allow indels/frameshifts with premature stop codons into clean cohort.") parser.add_argument("--full-length-only", action="store_true", help="Only include accessions that cover 100%% of query length.") parser.add_argument("--max-records", type=int, default=None, help="Stop after processing N lines (useful for testing/previews).") parser.add_argument("--sample-cap", type=int, default=10, help="Maximum sample accessions to record per haplotype in TSV (default: 10).") parser.add_argument("-o", "--out-prefix", "--out-dir", dest="out_prefix", type=str, default=None, help="Output prefix or directory for generated files (default: msa/).") parser.add_argument("--no-a2m", action="store_true", help="Disable generating .a2m alignment files.") return parser.parse_args() def main(): args = parse_args() if not args.dataset and not args.file: sys.exit("Error: Must specify either -d/--dataset or -f/-i/--file (or '-' for stdin).") stream_lexicmap_to_msa( ref_path=args.ref, dataset_id=args.dataset, file_path=args.file, base_url=args.base_url, api_key=args.api_key, api_key_file=args.api_key_file, min_coverage=args.min_coverage, min_coverage_partial=args.min_coverage_partial, min_pident=args.min_pident, max_internal_stops=args.max_internal_stops, allow_frameshifts=args.allow_frameshifts, full_length_only=args.full_length_only, max_records=args.max_records, sample_cap=args.sample_cap, out_prefix=args.out_prefix, write_a2m=not args.no_a2m ) if __name__ == "__main__": main() shell_command: >- python lexicmap_streamer.py -q '$(inputs.reference_fasta.path)' -f '$(inputs.lexicmap_results.path)' -o result --min-coverage $(inputs.min_coverage) --min-coverage-partial $(inputs.min_coverage_partial) --min-pident $(inputs.min_pident) --max-internal-stops $(inputs.max_internal_stops) --sample-cap $(inputs.sample_cap) --no-a2m $(inputs.allow_frameshifts ? '--allow-frameshifts' : '') outputs: - name: clean_msa_fasta type: data format: fasta from_work_dir: result.clean.msa.fasta label: Clean full-length haplotypes (MSA FASTA) - name: clean_haplotypes_tsv type: data format: tabular from_work_dir: result.clean.haplotypes.tsv label: Clean haplotype catalog (mutations, frequencies, samples) - name: clean_accessions_fasta type: data format: fasta from_work_dir: result.clean.accessions.fasta label: Clean full-length accessions (uncollapsed) - name: clean_expanded_accessions_fasta type: data format: fasta from_work_dir: result.clean_expanded.accessions.fasta label: Clean expanded accessions (full-length + partial) - name: flagged_accessions_fasta type: data format: fasta from_work_dir: result.flagged.accessions.fasta label: Flagged accessions (premature stops or low coverage) - name: flagged_tsv type: data format: tabular from_work_dir: result.flagged.tsv label: Flagged accession audit table - name: cohort_ledger_tsv type: data format: tabular from_work_dir: result.cohort_ledger.tsv label: Complete per-accession cohort ledger - name: summary_json type: data format: json from_work_dir: result.summary.json label: Per-gene tiling/QC/cohort summary (quantiles, counts, entropy)