diff --git a/README.md b/README.md index d4c5c5e..b53b940 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,7 @@ The repository contains a neutral simulated demo only. Raw reads are processed l - Optional YAML/UI `donor_ID`, same-donor contamination protection, and a donor-level combined collapsed/functional alignment and phylogeny workbench with sample-colored tips. - A bundled Alivibe pop-out editor with permissive biological-edit warnings, validated return, explicit tree recalculation, separately persisted alignment/frame/tree edits, and a detailed append-only edit audit. - Interactive UMI, artefact, agreement, MDS/APOBEC, and dinucleotide figures with labelled axes and SVG export. -- Reference-clipped codon-aware nucleotide/protein alignments for every functional-filtered collapsed variant. Sample proteins are aligned first; the reference is then added to that fixed profile without rearranging sample rows and is retained first under its supplied FASTA name. Passing variants are labelled `sample_vN_K rm=0.yy`, retaining the family-count suffix. Protein rows are verified before backtranslation; a sequence-preserving fallback repairs a structurally invalid MSA, while an unrecoverable functional-filter error is isolated to that sample and does not stop other samples. Sequences without a complete start-to-stop ORF are rejected. +- Codon-aware nucleotide/protein alignments for every functional-filtered collapsed variant. Functional screening aligns eligible sample ORFs, adds the reference to that fixed profile, and clips to the reference endpoints. Passing trimmed proteins are then realigned from scratch with bundled MAFFT 7.520 FFT-NS-2; MAFFT adds the reference to the passing-only profile without rearranging sample columns. Export retains every residue of these trimmed proteins, including terminal residues outside the newly aligned reference endpoints, and places the reference first under its supplied FASTA name. Passing variants retain `sample_vN_K rm=0.yy`; `rm` is the screening match score, before export realignment. Sequence preservation is checked before codon backtranslation, and a functional-filter error remains isolated to its sample. Sequences without a complete start-to-stop ORF are rejected. - A single compressed `.webporpid` result file, component exports, and a complete gzip-compressed tar bundle. - A subtle package-derived version label in the page header, so a deployed build can be identified immediately. diff --git a/cli-src/direct-mafft.mjs b/cli-src/direct-mafft.mjs index 1502e50..a77f119 100644 --- a/cli-src/direct-mafft.mjs +++ b/cli-src/direct-mafft.mjs @@ -2,40 +2,13 @@ import { readFile } from "node:fs/promises"; import { pathToFileURL } from "node:url"; import { Worker as NodeWorker } from "node:worker_threads"; -const FFT_NS_2_ARGUMENTS = [ - "-q", "0", "-E", "2", "-V", "-1.53", "-s", "0.0", "-W", "6", "-O", - "-C", "0-0", "-D", "-b", "62", "-g", "0", "-f", "-1.53", "-Q", "100.0", - "-h", "0", "-F", "-X", "0.1", "-x", "-1", "-i", "/input.fa", -]; - -const numericFasta = (sequences) => sequences.map((sequence, index) => `>${index}\n${sequence}\n`).join(""); - -function parseNumericAlignment(source, expected) { - const rows = new Map(); let index = -1, sequence = ""; - const finish = () => { if (index >= 0) { if (rows.has(index)) throw new Error("MAFFT returned a duplicate sequence identifier."); rows.set(index, sequence.toUpperCase()); } }; - for (const raw of source.split(/\r?\n/)) { - const line = raw.trim(); if (!line) continue; - if (line.startsWith(">")) { - finish(); index = Number(line.slice(1).trim()); sequence = ""; - if (!Number.isSafeInteger(index) || index < 0 || index >= expected.length) throw new Error("MAFFT returned an unknown sequence identifier."); - } else { if (index < 0) throw new Error("MAFFT returned sequence data before its first header."); sequence += line.replace(/\s/g, ""); } - } - finish(); - if (rows.size !== expected.length) throw new Error("MAFFT returned the wrong number of sequences."); - const aligned = expected.map((input, row) => { - const output = rows.get(row); if (output.replaceAll("-", "") !== input.toUpperCase()) throw new Error(`MAFFT changed sequence ${row + 1}.`); return output; - }); - const width = aligned[0]?.length ?? 0; - if (!width || aligned.some((row) => row.length !== width)) throw new Error("MAFFT returned a non-rectangular alignment."); - return aligned; -} +import { mafftArguments, numericMafftFasta, parseMafftAlignment, mafftProfileAddition } from "../src/mafft-msa-codec.ts"; export function createDirectMafftRunner(javascriptPath, wasmPath) { const factoryPromise = import(pathToFileURL(javascriptPath).href).then((module) => module.default); const wasmPromise = readFile(wasmPath); - return async (sequences, signal, _iterations = 0, scoringMode = "nucleotide", onProgress) => { + const run = async (sequences, signal, _iterations = 0, scoringMode = "nucleotide", onProgress, profileRows) => { if (signal?.aborted) throw new Error("Analysis cancelled."); - if (scoringMode === "amino-acid") throw new Error("The reference-panel MAFFT runner expects nucleotide sequences."); if (sequences.length < 2) return [...sequences]; const input = sequences.map((sequence) => String(sequence).toUpperCase()), stdout = []; if (input.some((sequence) => !sequence.length || /[^A-Z?*.-]/.test(sequence))) throw new Error("MAFFT input contains an unsupported symbol."); @@ -47,11 +20,14 @@ export function createDirectMafftRunner(javascriptPath, wasmPath) { lastProgress = now; onProgress?.({ detail }); } } }); - runtime.FS.writeFile("/input.fa", new TextEncoder().encode(numericFasta(input))); - try { const status = runtime.callMain([...FFT_NS_2_ARGUMENTS]); if (status) throw new Error(`MAFFT exited with status ${status}.`); } + runtime.FS.writeFile("/input.fa", new TextEncoder().encode(numericMafftFasta(input, scoringMode))); + try { const status = runtime.callMain(mafftArguments(scoringMode, profileRows)); if (status) throw new Error(`MAFFT exited with status ${status}.`); } finally { try { runtime.FS.unlink("/input.fa"); } catch { /* best effort */ } } - return parseNumericAlignment(stdout.join("\n"), input); + return parseMafftAlignment(stdout.join("\n"), input, scoringMode); }; + run.addToProfile = async (profileRows, sequence, signal) => mafftProfileAddition(profileRows, + await run([...profileRows, sequence], signal, 0, "amino-acid", undefined, profileRows.length)); + return run; } class MafftWorkerClient { @@ -71,7 +47,7 @@ class MafftWorkerClient { close() { return this.worker.terminate(); } } -/** Run each sample's batch MAFFT alignment on an isolated CPU worker. */ +/** Run each sample's MAFFT alignment on an isolated CPU worker. */ export function createMafftRunner(javascriptPath, wasmPath, size = 1, workerPath = new URL("../porpid-mafft-worker.mjs", import.meta.url)) { const count = Math.max(1, Math.floor(size)); const clients = Array.from({ length: count }, () => new MafftWorkerClient(new NodeWorker(workerPath))); let cursor = 0; @@ -79,6 +55,10 @@ export function createMafftRunner(javascriptPath, wasmPath, size = 1, workerPath if (signal?.aborted) throw new Error("Analysis cancelled."); return clients[cursor++ % clients.length].call({ javascriptPath, wasmPath, sequences: sequences.map(String), iterations, scoringMode }, onProgress); }; + run.addToProfile = async (profileRows, sequence, signal) => { + if (signal?.aborted) throw new Error("Analysis cancelled."); + return clients[cursor++ % clients.length].call({ javascriptPath, wasmPath, profileRows: profileRows.map(String), sequence }); + }; run.close = async () => { await Promise.all(clients.map((client) => client.close())); }; return run; } diff --git a/cli-src/porpid-cli.mjs b/cli-src/porpid-cli.mjs index 516b91f..67dee78 100644 --- a/cli-src/porpid-cli.mjs +++ b/cli-src/porpid-cli.mjs @@ -318,7 +318,8 @@ async function runPipeline({ inputPath, configPath, outputPath, workers, assets, const panelFilterRunner = createIndependentPanelFilterRunner(assets.panelWorkerPath); const downstreamStarted = performance.now(); let downstream; try { downstream = await postprocess(consensuses, contamination, config, undefined, msaRunner, workers, undefined, - { panelMsa: panelMsaRunner, panelFilter: panelFilterRunner }); } + { panelMsa: panelMsaRunner, panelFilter: panelFilterRunner, + functionalAlignment: { msa: panelMsaRunner, addToProfile: panelMsaRunner.addToProfile } }); } finally { await Promise.all([msaRunner.close?.(), panelMsaRunner.close?.()]); } downstream.summaries.forEach((summary, index) => { summary.demultiplexedReads = quality.perSample[index] ?? 0; diff --git a/cli-src/porpid-mafft-worker.mjs b/cli-src/porpid-mafft-worker.mjs index e06ec71..86e92de 100644 --- a/cli-src/porpid-mafft-worker.mjs +++ b/cli-src/porpid-mafft-worker.mjs @@ -6,7 +6,7 @@ parentPort.on("message", async (message) => { try { const nextKey = `${message.javascriptPath}\0${message.wasmPath}`; if (!runner || key !== nextKey) { key = nextKey; runner = createDirectMafftRunner(message.javascriptPath, message.wasmPath); } - const result = await runner(message.sequences, undefined, message.iterations, message.scoringMode, + const result = message.profileRows ? await runner.addToProfile(message.profileRows, message.sequence) : await runner(message.sequences, undefined, message.iterations, message.scoringMode, ({ detail }) => parentPort.postMessage({ id: message.id, progress: detail })); parentPort.postMessage({ id: message.id, result }); } catch (cause) { parentPort.postMessage({ id: message.id, error: cause instanceof Error ? cause.message : String(cause) }); } diff --git a/cli/chunks/direct-mafft-W7tDpHgK.mjs b/cli/chunks/direct-mafft-BgTc0bqp.mjs similarity index 58% rename from cli/chunks/direct-mafft-W7tDpHgK.mjs rename to cli/chunks/direct-mafft-BgTc0bqp.mjs index bc1a14d..e3f1de1 100644 --- a/cli/chunks/direct-mafft-W7tDpHgK.mjs +++ b/cli/chunks/direct-mafft-BgTc0bqp.mjs @@ -1,49 +1,60 @@ import { readFile } from "node:fs/promises"; import { Worker } from "node:worker_threads"; import { pathToFileURL } from "node:url"; -//#region cli-src/direct-mafft.mjs -const FFT_NS_2_ARGUMENTS = [ - "-q", - "0", - "-E", - "2", - "-V", - "-1.53", - "-s", - "0.0", - "-W", - "6", - "-O", - "-C", - "0-0", - "-D", - "-b", - "62", - "-g", - "0", - "-f", - "-1.53", - "-Q", - "100.0", - "-h", - "0", - "-F", - "-X", - "0.1", - "-x", - "-1", - "-i", - "/input.fa" -]; -const numericFasta = (sequences) => sequences.map((sequence, index) => `>${index}\n${sequence}\n`).join(""); -function parseNumericAlignment(source, expected) { +//#region src/mafft-msa-codec.ts +/** disttbfast flags used by MAFFT 7.520's FFT-NS-2 and --add wrappers. */ +function mafftArguments(mode, profileRows) { + return [ + "-q", + "0", + "-E", + "2", + "-V", + "-1.53", + "-s", + "0.0", + "-W", + "6", + "-O", + "-C", + "0-0", + mode === "amino-acid" ? "-P" : "-D", + "-b", + "62", + "-g", + "0", + "-f", + "-1.53", + "-Q", + "100.0", + "-h", + "0", + "-F", + "-X", + "0.1", + "-x", + "-1", + ...profileRows === void 0 ? [] : [ + "-K", + "-I", + "1" + ], + "-i", + "/input.fa" + ]; +} +function numericMafftFasta(sequences, mode) { + return sequences.map((sequence, index) => { + return `>${index}\n${mode === "amino-acid" ? sequence.replaceAll("*", "X") : sequence}\n`; + }).join(""); +} +function parseMafftAlignment(source, expected, mode) { const rows = /* @__PURE__ */ new Map(); let index = -1, sequence = ""; const finish = () => { - if (index >= 0) { - if (rows.has(index)) throw new Error("MAFFT returned a duplicate sequence identifier."); - rows.set(index, sequence.toUpperCase()); - } + if (index < 0) return; + if (rows.has(index)) throw new Error("MAFFT returned a duplicate sequence identifier."); + rows.set(index, sequence.toUpperCase()); }; for (const raw of source.split(/\r?\n/)) { const line = raw.trim(); @@ -61,20 +72,37 @@ function parseNumericAlignment(source, expected) { finish(); if (rows.size !== expected.length) throw new Error("MAFFT returned the wrong number of sequences."); const aligned = expected.map((input, row) => { - const output = rows.get(row); - if (output.replaceAll("-", "") !== input.toUpperCase()) throw new Error(`MAFFT changed sequence ${row + 1}.`); - return output; + const original = input.replaceAll("-", "").toUpperCase(), output = rows.get(row); + const normalized = mode === "amino-acid" ? original.replaceAll("*", "X") : original; + if (output.replaceAll("-", "") !== normalized) throw new Error(`MAFFT changed sequence ${row + 1}.`); + let cursor = 0; + return [...output].map((residue) => residue === "-" ? "-" : original[cursor++]).join(""); }); const width = aligned[0]?.length ?? 0; if (!width || aligned.some((row) => row.length !== width)) throw new Error("MAFFT returned a non-rectangular alignment."); return aligned; } +function mafftProfileAddition(profile, aligned) { + if (aligned.length !== profile.length + 1) throw new Error("MAFFT profile addition returned the wrong number of rows."); + const sampleRows = aligned.slice(0, profile.length); + const withoutGapColumns = (rows) => { + const columns = Array.from({ length: rows[0]?.length ?? 0 }, (_, column) => column).filter((column) => rows.some((row) => row[column] !== "-")); + return rows.map((row) => columns.map((column) => row[column]).join("")); + }; + const before = withoutGapColumns(profile), after = withoutGapColumns(sampleRows); + if (before.some((row, index) => row !== after[index])) throw new Error("MAFFT reference addition changed the existing sample protein alignment."); + return { + profileRows: sampleRows, + sequence: aligned[profile.length] + }; +} +//#endregion +//#region cli-src/direct-mafft.mjs function createDirectMafftRunner(javascriptPath, wasmPath) { const factoryPromise = import(pathToFileURL(javascriptPath).href).then((module) => module.default); const wasmPromise = readFile(wasmPath); - return async (sequences, signal, _iterations = 0, scoringMode = "nucleotide", onProgress) => { + const run = async (sequences, signal, _iterations = 0, scoringMode = "nucleotide", onProgress, profileRows) => { if (signal?.aborted) throw new Error("Analysis cancelled."); - if (scoringMode === "amino-acid") throw new Error("The reference-panel MAFFT runner expects nucleotide sequences."); if (sequences.length < 2) return [...sequences]; const input = sequences.map((sequence) => String(sequence).toUpperCase()), stdout = []; if (input.some((sequence) => !sequence.length || /[^A-Z?*.-]/.test(sequence))) throw new Error("MAFFT input contains an unsupported symbol."); @@ -92,17 +120,19 @@ function createDirectMafftRunner(javascriptPath, wasmPath) { } } }); - runtime.FS.writeFile("/input.fa", new TextEncoder().encode(numericFasta(input))); + runtime.FS.writeFile("/input.fa", new TextEncoder().encode(numericMafftFasta(input, scoringMode))); try { - const status = runtime.callMain([...FFT_NS_2_ARGUMENTS]); + const status = runtime.callMain(mafftArguments(scoringMode, profileRows)); if (status) throw new Error(`MAFFT exited with status ${status}.`); } finally { try { runtime.FS.unlink("/input.fa"); } catch {} } - return parseNumericAlignment(stdout.join("\n"), input); + return parseMafftAlignment(stdout.join("\n"), input, scoringMode); }; + run.addToProfile = async (profileRows, sequence, signal) => mafftProfileAddition(profileRows, await run([...profileRows, sequence], signal, 0, "amino-acid", void 0, profileRows.length)); + return run; } var MafftWorkerClient = class { constructor(worker) { @@ -146,7 +176,7 @@ var MafftWorkerClient = class { return this.worker.terminate(); } }; -/** Run each sample's batch MAFFT alignment on an isolated CPU worker. */ +/** Run each sample's MAFFT alignment on an isolated CPU worker. */ function createMafftRunner(javascriptPath, wasmPath, size = 1, workerPath = new URL("../porpid-mafft-worker.mjs", import.meta.url)) { const count = Math.max(1, Math.floor(size)); const clients = Array.from({ length: count }, () => new MafftWorkerClient(new Worker(workerPath))); @@ -161,10 +191,19 @@ function createMafftRunner(javascriptPath, wasmPath, size = 1, workerPath = new scoringMode }, onProgress); }; + run.addToProfile = async (profileRows, sequence, signal) => { + if (signal?.aborted) throw new Error("Analysis cancelled."); + return clients[cursor++ % clients.length].call({ + javascriptPath, + wasmPath, + profileRows: profileRows.map(String), + sequence + }); + }; run.close = async () => { await Promise.all(clients.map((client) => client.close())); }; return run; } //#endregion -export { createMafftRunner as n, createDirectMafftRunner as t }; +export { createMafftRunner as n, mafftProfileAddition as r, createDirectMafftRunner as t }; diff --git a/cli/chunks/direct-msa-Bnko_59s.mjs b/cli/chunks/direct-msa-gemMlFDe.mjs similarity index 99% rename from cli/chunks/direct-msa-Bnko_59s.mjs rename to cli/chunks/direct-msa-gemMlFDe.mjs index ae91a75..44b1ca6 100644 --- a/cli/chunks/direct-msa-Bnko_59s.mjs +++ b/cli/chunks/direct-msa-gemMlFDe.mjs @@ -1,4 +1,4 @@ -import { t as WASI } from "./dist-BPeWFyil.mjs"; +import { t as WASI } from "./dist-DaUuEvmK.mjs"; import { readFile } from "node:fs/promises"; import { Worker } from "node:worker_threads"; //#region src/alivibe-msa-codec.ts diff --git a/cli/chunks/dist-BPeWFyil.mjs b/cli/chunks/dist-DaUuEvmK.mjs similarity index 100% rename from cli/chunks/dist-BPeWFyil.mjs rename to cli/chunks/dist-DaUuEvmK.mjs diff --git a/cli/chunks/independent-panel-filter-DTwAa3St.mjs b/cli/chunks/independent-panel-filter-DxKQKnbW.mjs similarity index 100% rename from cli/chunks/independent-panel-filter-DTwAa3St.mjs rename to cli/chunks/independent-panel-filter-DxKQKnbW.mjs diff --git a/cli/chunks/wasm-runtime-D2YnsFQQ.mjs b/cli/chunks/wasm-runtime-Che1h7b0.mjs similarity index 99% rename from cli/chunks/wasm-runtime-D2YnsFQQ.mjs rename to cli/chunks/wasm-runtime-Che1h7b0.mjs index bbf64a2..f5f5037 100644 --- a/cli/chunks/wasm-runtime-D2YnsFQQ.mjs +++ b/cli/chunks/wasm-runtime-Che1h7b0.mjs @@ -1,5 +1,5 @@ import { n as BinaryWriter, t as BinaryReader } from "./binary-68r8u1WT.mjs"; -import { t as WASI } from "./dist-BPeWFyil.mjs"; +import { t as WASI } from "./dist-DaUuEvmK.mjs"; //#region src/wasm-runtime.ts const encoder = new TextEncoder(); const decoder = new TextDecoder(); diff --git a/cli/porpid-cli.mjs b/cli/porpid-cli.mjs index fc01776..e258aea 100755 --- a/cli/porpid-cli.mjs +++ b/cli/porpid-cli.mjs @@ -1,9 +1,9 @@ #!/usr/bin/env node import { a as resolveReferenceFiles, i as parseConfigYaml, n as createFastTreeRunner, o as resultConfig, r as compileConfig } from "./chunks/direct-fasttree-CuMSsXhU.mjs"; -import { a as encodeAlivibeMsaSequences, i as decodeAlivibeMsaSequences, n as createMsaRunner, r as assertAlivibeMsaResult } from "./chunks/direct-msa-Bnko_59s.mjs"; -import { n as filterQueriesAgainstPanel, t as addSequenceToProfile } from "./chunks/independent-panel-filter-DTwAa3St.mjs"; -import { a as encodeFamilyModel, c as mergeFamilyCounts, i as decodeFamilyModel, l as mergeStats, n as decodeConsensusOutput, o as makeCutoffValues, r as decodeFamilyCounts, s as makeCutoffs } from "./chunks/wasm-runtime-D2YnsFQQ.mjs"; -import { n as createMafftRunner } from "./chunks/direct-mafft-W7tDpHgK.mjs"; +import { a as encodeAlivibeMsaSequences, i as decodeAlivibeMsaSequences, n as createMsaRunner, r as assertAlivibeMsaResult } from "./chunks/direct-msa-gemMlFDe.mjs"; +import { n as createMafftRunner, r as mafftProfileAddition } from "./chunks/direct-mafft-BgTc0bqp.mjs"; +import { n as filterQueriesAgainstPanel, t as addSequenceToProfile } from "./chunks/independent-panel-filter-DxKQKnbW.mjs"; +import { a as encodeFamilyModel, c as mergeFamilyCounts, i as decodeFamilyModel, l as mergeStats, n as decodeConsensusOutput, o as makeCutoffValues, r as decodeFamilyCounts, s as makeCutoffs } from "./chunks/wasm-runtime-Che1h7b0.mjs"; import { createHash } from "node:crypto"; import { createReadStream } from "node:fs"; import { mkdir, mkdtemp, open, readFile, rm, stat, writeFile } from "node:fs/promises"; @@ -647,8 +647,22 @@ async function runIndependentPanelFilter(sequences, panelRows, signal, workers = //#region src/mafft-msa-runtime.ts /** Run MAFFT 7.520 FFT-NS-2 in an isolated worker with growable WASM memory. */ async function runMafftFftnsMsa(sequences, signal, _iterations = 0, scoringMode = "nucleotide", onProgress) { + if (signal?.aborted) throw new DOMException("MAFFT alignment cancelled.", "AbortError"); if (sequences.length < 2) return [...sequences]; - if (scoringMode === "amino-acid") throw new Error("The reference-panel MAFFT runner expects nucleotide sequences."); + return runMafft({ + sequences: sequences.map(String), + scoringMode + }, signal, onProgress); +} +/** Equivalent to mafft --amino --add reference.fa sample-alignment.fa. */ +async function runMafftProfileAdd(profileRows, sequence, signal) { + return mafftProfileAddition(profileRows, await runMafft({ + sequences: [...profileRows, sequence], + scoringMode: "amino-acid", + profileRows: profileRows.length + }, signal)); +} +async function runMafft(request, signal, onProgress) { const worker = new Worker(new URL("./mafft-msa-worker.ts", import.meta.url), { type: "module" }); return new Promise((resolve, reject) => { let settled = false; @@ -660,11 +674,11 @@ async function runMafftFftnsMsa(sequences, signal, _iterations = 0, scoringMode return true; }; const abort = () => { - if (finish()) reject(new DOMException("MAFFT panel alignment cancelled.", "AbortError")); + if (finish()) reject(new DOMException("MAFFT alignment cancelled.", "AbortError")); }; worker.onmessage = (event) => { if (event.data.type === "progress") { - onProgress?.({ detail: event.data.detail ?? "MAFFT is aligning candidate sequences" }); + onProgress?.({ detail: event.data.detail ?? "MAFFT is aligning sequences" }); return; } if (!finish()) return; @@ -679,7 +693,7 @@ async function runMafftFftnsMsa(sequences, signal, _iterations = 0, scoringMode return; } signal?.addEventListener("abort", abort, { once: true }); - worker.postMessage({ sequences: sequences.map(String) }); + worker.postMessage(request); }); } //#endregion @@ -1199,6 +1213,10 @@ async function runScalableMsa(sequences, runMsa, signal, iterations = 3, scoring } //#endregion //#region src/postprocess.ts +const defaultFunctionalAlignment = { + msa: runMafftFftnsMsa, + addToProfile: runMafftProfileAdd +}; const degap = (sequence) => sequence.replaceAll("-", "").toUpperCase(); const fasta$1 = (rows) => rows.map((row) => `>${row.name}\n${row.sequence.match(/.{1,80}/g)?.join("\n") ?? ""}`).join("\n") + (rows.length ? "\n" : ""); function quantile(values, probability) { @@ -1515,7 +1533,7 @@ function preserveProteinSequences(inputs, aligned) { * Projecting every resulting gap back as a codon triplet prevents an aligner * from manufacturing one- or two-base frameshifts. */ -async function functionalFilterBatch(referenceName, reference, sequences, threshold, runMsa, signal) { +async function functionalFilterBatch(referenceName, reference, sequences, threshold, runMsa, signal, exportAlignment = defaultFunctionalAlignment) { const outcomes = Array(sequences.length), coding = []; sequences.forEach((raw, index) => { const sequence = degap(raw); @@ -1570,12 +1588,28 @@ async function functionalFilterBatch(referenceName, reference, sequences, thresh referenceMatch: Number(rawRatio.toFixed(2)) }; }); - return { + const passing = outcomes.filter((outcome) => outcome.passed); + if (!passing.length) return { outcomes, referenceName: resolvedReferenceName, referenceNt: referenceRegion, referenceAa: translateAlignedNucleotides(referenceRegion, 0) }; + const proteins = passing.map((outcome) => outcome.aa); + const profile = proteins.length > 1 ? await runScalableMsa(proteins, exportAlignment.msa, signal, 0, "amino-acid") : proteins; + if (profile.length !== proteins.length || profile.some((row, index) => row.length !== profile[0].length || degap(row) !== proteins[index])) throw new Error("The functional export alignment did not preserve its input proteins."); + const aligned = await exportAlignment.addToProfile(profile, translate(referenceCoding), signal); + if (aligned.profileRows.length !== passing.length || degap(aligned.sequence) !== translate(referenceCoding) || aligned.profileRows.some((row, index) => row.length !== aligned.sequence.length || degap(row) !== proteins[index])) throw new Error("The functional export reference addition did not preserve its input proteins."); + passing.forEach((outcome, index) => { + outcome.alignedAa = aligned.profileRows[index]; + outcome.alignedNt = backtranslate(outcome.alignedAa, outcome.nt); + }); + return { + outcomes, + referenceName: resolvedReferenceName, + referenceNt: backtranslate(aligned.sequence, referenceCoding), + referenceAa: aligned.sequence + }; } async function postprocess(consensuses, contamination, config, signal, runMsa = runAlivibeMsa, sampleConcurrency = 1, onProgress, options = {}) { const panelMsa = options.panelMsa ?? runMafftFftnsMsa; @@ -1734,10 +1768,10 @@ async function postprocess(consensuses, contamination, config, signal, runMsa = functionalFilterErrors: Object.assign({}, ...completed.map((output) => output.functionalFilterErrors)), collapseSeconds: 0 }; - return options.collapse === false ? combined : collapsePostprocess(combined, config, signal, options.onCollapseProgress, runMsa); + return options.collapse === false ? combined : collapsePostprocess(combined, config, signal, options.onCollapseProgress, runMsa, void 0, options.functionalAlignment); } /** Run family-count-preserving haplotype collapse as its own resumable stage. */ -async function collapsePostprocess(output, config, signal, onProgress, runMsa = runAlivibeMsa, sampleNames) { +async function collapsePostprocess(output, config, signal, onProgress, runMsa = runAlivibeMsa, sampleNames, functionalAlignment = defaultFunctionalAlignment) { const started = performance.now(), alignments = { ...output.alignments }, referenceAlignments = { ...output.referenceAlignments }; const collapseGroups = { ...output.collapseGroups }, summaries = output.summaries.map((summary) => ({ ...summary })); const functionalFilterErrors = { ...output.functionalFilterErrors }; @@ -1775,7 +1809,7 @@ async function collapsePostprocess(output, config, signal, onProgress, runMsa = }); let batch; try { - batch = await functionalFilterBatch(sample.functionalReferenceSequence.name, sample.functionalReferenceSequence.sequence, collapsedRows.map((row) => row.sequence), sample.functionalMatchOverride ?? config.parameters.functionalMatchThreshold, runMsa, signal); + batch = await functionalFilterBatch(sample.functionalReferenceSequence.name, sample.functionalReferenceSequence.sequence, collapsedRows.map((row) => row.sequence), sample.functionalMatchOverride ?? config.parameters.functionalMatchThreshold, runMsa, signal, functionalAlignment); } catch (cause) { if (signal?.aborted || cause instanceof DOMException && cause.name === "AbortError") throw cause; const message = cause instanceof Error ? cause.message : String(cause); @@ -5179,7 +5213,11 @@ async function runPipeline({ inputPath, configPath, outputPath, workers, assets, try { downstream = await postprocess(consensuses, contamination, config, void 0, msaRunner, workers, void 0, { panelMsa: panelMsaRunner, - panelFilter: panelFilterRunner + panelFilter: panelFilterRunner, + functionalAlignment: { + msa: panelMsaRunner, + addToProfile: panelMsaRunner.addToProfile + } }); } finally { await Promise.all([msaRunner.close?.(), panelMsaRunner.close?.()]); diff --git a/cli/porpid-mafft-worker.mjs b/cli/porpid-mafft-worker.mjs index b45aa2b..25c6fe3 100644 --- a/cli/porpid-mafft-worker.mjs +++ b/cli/porpid-mafft-worker.mjs @@ -1,4 +1,4 @@ -import { t as createDirectMafftRunner } from "./chunks/direct-mafft-W7tDpHgK.mjs"; +import { t as createDirectMafftRunner } from "./chunks/direct-mafft-BgTc0bqp.mjs"; import { parentPort } from "node:worker_threads"; //#region cli-src/porpid-mafft-worker.mjs let key = "", runner; @@ -9,7 +9,7 @@ parentPort.on("message", async (message) => { key = nextKey; runner = createDirectMafftRunner(message.javascriptPath, message.wasmPath); } - const result = await runner(message.sequences, void 0, message.iterations, message.scoringMode, ({ detail }) => parentPort.postMessage({ + const result = message.profileRows ? await runner.addToProfile(message.profileRows, message.sequence) : await runner(message.sequences, void 0, message.iterations, message.scoringMode, ({ detail }) => parentPort.postMessage({ id: message.id, progress: detail })); diff --git a/cli/porpid-msa-worker.mjs b/cli/porpid-msa-worker.mjs index 50e0db7..c58542d 100644 --- a/cli/porpid-msa-worker.mjs +++ b/cli/porpid-msa-worker.mjs @@ -1,4 +1,4 @@ -import { t as createDirectMsaRunner } from "./chunks/direct-msa-Bnko_59s.mjs"; +import { t as createDirectMsaRunner } from "./chunks/direct-msa-gemMlFDe.mjs"; import { parentPort } from "node:worker_threads"; //#region cli-src/porpid-msa-worker.mjs let wasmPath, runner; diff --git a/cli/porpid-panel-worker.mjs b/cli/porpid-panel-worker.mjs index d311b89..05af67f 100644 --- a/cli/porpid-panel-worker.mjs +++ b/cli/porpid-panel-worker.mjs @@ -1,4 +1,4 @@ -import { n as filterQueriesAgainstPanel } from "./chunks/independent-panel-filter-DTwAa3St.mjs"; +import { n as filterQueriesAgainstPanel } from "./chunks/independent-panel-filter-DxKQKnbW.mjs"; import { parentPort } from "node:worker_threads"; //#region cli-src/porpid-panel-worker.mjs parentPort.on("message", (message) => { diff --git a/cli/porpid-worker.mjs b/cli/porpid-worker.mjs index 37c689d..98d2ed8 100644 --- a/cli/porpid-worker.mjs +++ b/cli/porpid-worker.mjs @@ -1,4 +1,4 @@ -import { t as WebPorpidRuntime } from "./chunks/wasm-runtime-D2YnsFQQ.mjs"; +import { t as WebPorpidRuntime } from "./chunks/wasm-runtime-Che1h7b0.mjs"; import { readFile } from "node:fs/promises"; import { parentPort } from "node:worker_threads"; //#region cli-src/porpid-worker.mjs diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index de91ad0..6a967a0 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -56,7 +56,7 @@ MSA is monolithic up to 8,000 sequences and 128 MiB of input residues. Beyond ei After non-functional filtering, identical ungapped haplotypes are collapsed. Variants are ranked by decreasing member-family count (lexical haplotype tie-break) and named `sample_vN_K`; uncollapsed FASTA/display labels are `sample_UMI fs=K minag=0.yy`, while the internal join key remains `sample_UMI`. Each collapsed variant stores those stable member identifiers, so bubble area is proportional to family count rather than read count. Minimum agreement remains attached only to each UMI family and is not assigned to the collapsed haplotype. FastTree consumes this collapsed nucleotide alignment by default. The uncollapsed alignment is retained for optional on-demand family-level inference and minimum-agreement colouring. The browser uses the packaged double-precision BioWASM/Aioli FastTree runtime; the CLI invokes the same packaged FastTree WASM directly. Browser tree failure is logged and produces an explicit zero-branch star fallback so a completed analysis remains loadable. -The functional filter runs only on collapsed variants. It rejects sequences without a complete start-to-stop ORF, translates the remaining ORFs, and first aligns the sample variants to one another. Each returned row is checked against its exact protein before backtranslation. Symbol-only normalization is restored; a residue-count mismatch switches to an anchor/profile construction that proves every ungapped row equals its input. The translated reference is then affine-aligned to that frozen sample profile: only shared all-gap columns may be inserted into existing sample rows, so adding the reference cannot rearrange their MSA. Gaps are back-translated as codon triplets. The stored functional-pass alignment is sliced from the first through last non-gap reference nucleotide and retains the reference as row one under its supplied FASTA name. Passing variants are displayed as `sample_vN_K rm=0.yy`, retaining identity, family count, and two-decimal match in one label. If both alignment routes fail, the failure is represented as sample-scoped state; no partial functional calls are retained and collapse/other samples continue. The alignment is retained as an editable view and can receive its own on-demand FastTree tree. +The functional filter runs only on collapsed variants. It rejects sequences without a complete start-to-stop ORF, translates the remaining ORFs, and first aligns the sample variants to one another. Each returned row is checked against its exact protein before backtranslation. Symbol-only normalization is restored; a residue-count mismatch switches to an anchor/profile construction that proves every ungapped row equals its input. The translated reference is then affine-aligned to that frozen sample profile: only shared all-gap columns may be inserted into existing sample rows, so adding the reference cannot rearrange their MSA. Gaps are back-translated as codon triplets. The screening alignment is clipped to the first/last non-gap reference nucleotide to produce each variant’s trimmed coding sequence. After screening, only passing trimmed proteins are degapped and freshly aligned with bundled MAFFT 7.520 FFT-NS-2. The translated reference is added to this fixed passing-only profile using MAFFT `--amino --add reference.fa passing-protein-alignment.fa`, then moved to row one under its supplied FASTA name. Sample columns and every ungapped residue are verified; terminal stop symbols are preserved and gaps are backtranslated as codon triplets. This export alignment is not clipped a second time, so it retains all residues of the already-trimmed passing sequences and the reference. Filter decisions and `rm` annotations remain based on the screening alignment. Stored projects need collapse/functional recomputation to obtain the new export layout. Passing variants are displayed as `sample_vN_K rm=0.yy`, retaining identity, family count, and two-decimal match in one label. If both alignment routes fail, the failure is represented as sample-scoped state; no partial functional calls are retained and collapse/other samples continue. The alignment is retained as an editable view and can receive its own on-demand FastTree tree. The viewer identifies the most frequent observed aligned nucleotide sequence, weighting collapsed representatives by retained UMI-family multiplicity and grouping identical uncollapsed tips on the fly. It roots the displayed tree on the branch endpoint immediately above a deterministic tip carrying that modal sequence; the modal-tip edge is exactly zero and the complete original pendant length is assigned to the opposite root edge. A midpoint-root action finds the weighted tree diameter in linear time and splits only the edge containing its midpoint. Regression tests require the complete tip set and every pairwise patristic distance to remain unchanged after either operation. diff --git a/docs/METHODS.md b/docs/METHODS.md index 9ba0494..981826f 100644 --- a/docs/METHODS.md +++ b/docs/METHODS.md @@ -55,7 +55,7 @@ One intentional reproducibility difference remains: Julia resolves ambiguous IUP - The artefact cutoff is `ceil(quantile(non-contaminant family sizes, q) * artefactFraction)` with per-sample overrides. - Agreement and contamination decisions are retained separately for audit and export. - Reference-panel filtering offers two alignment strategies before applying the same extraction and maximum-subarray panel score. The default bundles MAFFT 7.520 FFT-NS-2 to jointly align candidate families, converts that alignment to a candidate profile, and affine-profile-aligns it to the supplied fixed panel. The optional `independent-query` mode builds the fixed panel profile once per worker chunk and uses a length-scaled adaptive-banded affine alignment for each candidate separately. It removes query-to-query MSA work, preserves input order, distributes independent candidates across workers, and expands to the full dynamic-programming matrix whenever traceback approaches a band edge. Both modes retain the same panel-column expected-match costs, gap-open/extend behavior, extraction coordinates, indel-ignoring probability, and strict panel-score cutoff. -- The functional filter is evaluated only on exact collapsed nucleotide variants; uncollapsed UMI-family consensuses never receive functional decisions. It selects the longest complete start-to-stop ORF across three frames and rejects a sequence as `noORF-reject` when none exists. Eligible sample ORFs are translated and amino-acid aligned first. Every aligned row must preserve its exact input residues: benign alphabet normalization is restored onto the returned gap path, while a changed residue count triggers a deterministic anchor/profile fallback that cannot delete or invent residues. The reference is then affine-aligned to that fixed sample profile, inserting only shared all-gap columns into existing rows, and the result is backtranslated with complete codon gaps. The supplied reference FASTA name is retained on the first row and the alignment is clipped to its first/last non-gap nucleotide. Passing variants use `sample_vN_K rm=0.yy`, retaining the family-count suffix. An unrecoverable functional-alignment error is stored against only that sample; its functional outputs are omitted while collapse, other samples, and later stages continue. +- The functional filter is evaluated only on exact collapsed nucleotide variants; uncollapsed UMI-family consensuses never receive functional decisions. It selects the longest complete start-to-stop ORF across three frames and rejects a sequence as `noORF-reject` when none exists. Eligible sample ORFs are translated and amino-acid aligned first. Every aligned row must preserve its exact input residues: benign alphabet normalization is restored onto the returned gap path, while a changed residue count triggers a deterministic anchor/profile fallback that cannot delete or invent residues. The reference is then affine-aligned to that fixed sample profile, inserting only shared all-gap columns into existing rows, and the result is backtranslated with complete codon gaps. The screening alignment is clipped to the first/last non-gap reference nucleotide to produce each variant’s trimmed coding sequence. After screening, only passing trimmed proteins are degapped and freshly aligned with bundled MAFFT 7.520 FFT-NS-2. The translated reference is added to this fixed passing-only profile using MAFFT `--amino --add reference.fa passing-protein-alignment.fa`, then moved to row one under its supplied FASTA name. Sample columns and every ungapped residue are verified; terminal stop symbols are preserved and gaps are backtranslated as codon triplets. This export alignment is not clipped a second time, so it retains all residues of the already-trimmed passing sequences and the reference. Filter decisions and `rm` annotations remain based on the screening alignment. Stored projects need collapse/functional recomputation to obtain the new export layout. Passing variants use `sample_vN_K rm=0.yy`, retaining the family-count suffix. An unrecoverable functional-alignment error is stored against only that sample; its functional outputs are omitted while collapse, other samples, and later stages continue. - APOBEC summaries use the local four-state mutation grid and posterior integration; they are reporting fields rather than a filtering gate. - Accepted nucleotide sequences are aligned with the packaged MSA and nucleotide trees use packaged FastTree. Protein exploration does not depend on the optional functional filter: it directly translates every active aligned nucleotide row, with `--- → -` and mixed-gap/ambiguous/incomplete codons → `X`, using a persisted frame offset of 0, 1, or 2. diff --git a/docs/PARITY.md b/docs/PARITY.md index ac5d95c..2a4c53e 100644 --- a/docs/PARITY.md +++ b/docs/PARITY.md @@ -64,3 +64,7 @@ bash scripts/parity-demo.sh - a gzipped CLI analysis, result inspection, FASTA export, and identical consensus with one versus two workers. The parity fixtures are neutral synthetic strings. The bundled demo is simulated and contains no original biological demo material. + +## Functional export realignment + +The parity counts above describe functional screening. Functional-pass exports now rebuild the MSA from only passing, reference-trimmed proteins with MAFFT 7.520 FFT-NS-2 and add the reference to that fixed profile. Gap layouts therefore differ from the historical exports; decisions, screening match scores, and ungapped trimmed coding sequences remain unchanged. A rerun of the public nanopore demo retains all 335 donor_1_REN passes and reduces its reference-first protein alignment from 898 to 873 columns while preserving every residue and terminal stop. The manually reported 869-column alignment has not been reproduced without its MAFFT version and exact command/options. diff --git a/docs/RESULT_FORMAT.md b/docs/RESULT_FORMAT.md index c7ded3b..10988fd 100644 --- a/docs/RESULT_FORMAT.md +++ b/docs/RESULT_FORMAT.md @@ -23,7 +23,7 @@ The browser and CLI share the same encoder/decoder. The current decoder caps com - **Consensus records:** stable `sample_UMI` identifiers, sample, UMI, family size, nucleotide sequence, minimum agreement, and low-agreement site details. FASTA/display labels append `fs=K minag=0.yy` without changing the stored join key. - **Contamination calls:** one decision per consensus family, with nearest non-self label/distance and primary-discard versus wider-suspect status. A primary decision takes precedence over the wider suspect pass. - **Post-processing records:** family consensus/aligned sequences, every family-level filter decision, panel score, and optional APOBEC posterior summary. New results do not attach functional decisions to these uncollapsed records. -- **Per-sample products:** collapsed and annotated uncollapsed nucleotide FASTA alignments, reference-first/reference-clipped codon-aware functional-pass nucleotide/protein alignments using the supplied reference FASTA name, directly translated protein views, optional Newick trees, and aligned reference rows used for coordinate display. +- **Per-sample products:** collapsed and annotated uncollapsed nucleotide FASTA alignments, reference-first codon-aware functional-pass nucleotide/protein alignments using the supplied reference FASTA name, freshly MAFFT-aligned from passing trimmed proteins, directly translated protein views, optional Newick trees, and aligned reference rows used for coordinate display. - **Collapse membership:** abundance-ranked `sample_vN_K` representative, stable member identifiers, a multiplicity that counts UMI families rather than reads, and collapsed-only functional decision/trimmed outputs. Functional calls may store a two-decimal `referenceMatch`; passing functional headers are `sample_vN_K rm=0.yy`, preserving the count suffix. Minimum agreement remains exclusively attached to the member UMI families. - **Functional-filter errors:** optional sample-to-message entries identify a sample whose functional alignment could not be recovered. That sample has no partial functional decisions or functional alignment, while its collapse output and all other samples remain valid. - **Input mapping/run options/stage state:** exact uploaded-file to YAML-slot assignments, every requested defer option, automatic versus external scratch selection, an explicit completed/deferred/skipped record for contamination, post-processing, collapse, and tree stages, and whether stored post-processing applied or bypassed contamination. Scratch directory handles and paths are never serialized. A consensus-only project is valid; missing post-processing records are not interpreted as passes. Contamination may be incomplete while later stages are complete, but those outputs must be marked bypassed and cannot contain contamination rejects. @@ -54,7 +54,7 @@ Validation requires unique sample, family, consensus, and post-processing identi | `uncollapsed-nucleotide-alignment` | Every retained UMI-family consensus in its active alignment with `fs` and two-decimal `minag` annotations | | `uncollapsed-protein-alignment` | Direct translation of the active uncollapsed alignment and saved frame | | `uncollapsed-newick` | Optional on-demand family-level tree | -| `functional-nucleotide-alignment` | Functional reference as row one, then every passing `sample_vN_K rm=0.yy` variant, clipped to the reference endpoints | +| `functional-nucleotide-alignment` | Functional reference as row one, then every passing `sample_vN_K rm=0.yy` variant; passing trimmed sequences are realigned without further clipping (`rm` remains the screening score) | | `functional-protein-alignment` | Direct translation of that reference-first functional nucleotide alignment and saved frame | | `functional-newick` | Optional on-demand functional-sequence tree | | `log` | Persistent plain-text run log | diff --git a/src/Methods.tsx b/src/Methods.tsx index ad798ad..0f9455d 100644 --- a/src/Methods.tsx +++ b/src/Methods.tsx @@ -139,10 +139,11 @@ function DownstreamPage() { -
+

This filter runs only after exact nucleotide collapse and is evaluated once per collapsed variant. It is never applied separately to the uncollapsed UMI-family consensuses. Sequences containing symbols outside A/C/G/T receive ambiguousSymbols-reject. The longest complete start-to-stop ORF across three frames is selected; a sequence with no such ORF receives noORF-reject and is not aligned.

Eligible collapsed ORFs are translated and aligned to one another as amino acids first. Before backtranslation, every row is checked against its exact input protein. Harmless alphabet normalization is restored onto the aligner’s gap path. If the aligner changed a row’s residue count, a deterministic anchor/profile fallback realigns every protein while proving that no residue was deleted or invented. That biological sample MSA is then frozen. The translated functional reference is added by affine sequence-to-profile alignment: query insertions add shared all-gap columns to the sample rows, but existing sample residues and their column relationships cannot move. Every amino-acid gap is projected back as exactly ---, so alignment cannot introduce a one- or two-base frameshift. The first returned row uses the actual header from the supplied functional-reference FASTA. If even the sequence-preserving fallback fails, the error and absence of functional calls are stored for that sample only; its collapse output, other samples, and later pipeline work continue.

-

The nucleotide alignment is sliced from the first through last non-gap nucleotide of the functional reference, which is retained as its first row. Query-only terminal sequence is excluded from exports and tree inference. Reference match is exact aligned nucleotide identity divided by the clipped alignment width; the existing three-significant-digit value controls functionalMatchThreshold, while the stored/displayed annotation is rounded to two decimal places. A passing variant keeps its complete collapsed identity and is labelled sample_vN_K rm=0.yy in functional FASTA/alignment outputs.

+

The screening nucleotide alignment is sliced from the first through last non-gap nucleotide of the functional reference. Query-only terminal sequence is excluded when constructing each variant’s trimmed coding sequence. Reference match is exact aligned nucleotide identity divided by the clipped screening width; the existing three-significant-digit value controls functionalMatchThreshold, while the stored/displayed annotation is rounded to two decimal places. A passing variant keeps its complete collapsed identity and is labelled sample_vN_K rm=0.yy in functional FASTA/alignment outputs. This annotation remains the screening score, rather than a score recalculated from the export layout.

+

After these decisions, all gaps are removed from the passing trimmed proteins and a fresh bundled MAFFT 7.520 FFT-NS-2 amino-acid alignment is built. Rejected ORFs cannot leave gap columns in this final sample profile. The ungapped translated reference is then added with MAFFT --amino --add reference.fa passing-protein-alignment.fa: the passing alignment is the existing profile and the reference is the new sequence. The reference is moved to row one for export; sample order and existing column relationships are checked. Stop symbols are restored after alignment, and all gaps are backtranslated as codon triplets. This final alignment is not clipped again: all residues of each already-trimmed passing protein and of the supplied reference are preserved, even where the reference has terminal gaps. Functional decisions and rm scores are unchanged by this realignment. Existing saved projects retain their stored alignments until collapse/functional filtering is recomputed. Alignment width can vary with MAFFT strategy/version; retaining terminal stop symbols can also differ from a standalone MAFFT run that drops them.

diff --git a/src/mafft-msa-codec.ts b/src/mafft-msa-codec.ts new file mode 100644 index 0000000..094e820 --- /dev/null +++ b/src/mafft-msa-codec.ts @@ -0,0 +1,70 @@ +import type { ProfileAddition } from "./independent-panel-filter.ts"; + +export type MafftScoringMode = "literal" | "nucleotide" | "amino-acid"; +export interface MafftRequest { sequences: string[]; scoringMode: MafftScoringMode; profileRows?: number } + +/** disttbfast flags used by MAFFT 7.520's FFT-NS-2 and --add wrappers. */ +export function mafftArguments(mode: MafftScoringMode, profileRows?: number): string[] { + return [ + "-q", "0", "-E", "2", "-V", "-1.53", "-s", "0.0", "-W", "6", "-O", + "-C", "0-0", mode === "amino-acid" ? "-P" : "-D", "-b", "62", "-g", "0", "-f", "-1.53", "-Q", "100.0", + "-h", "0", "-F", "-X", "0.1", "-x", "-1", + // Existing profile rows come first; exactly one new reference comes last. + ...(profileRows === undefined ? [] : ["-K", "-I", "1"]), "-i", "/input.fa", + ]; +} + +export function numericMafftFasta(sequences: readonly string[], mode: MafftScoringMode): string { + return sequences.map((sequence, index) => { + // MAFFT drops '*' in protein input. Align it as X, then restore the exact + // original residue at that position before codon backtranslation. + const input = mode === "amino-acid" ? sequence.replaceAll("*", "X") : sequence; + return `>${index}\n${input}\n`; + }).join(""); +} + +export function parseMafftAlignment(source: string, expected: readonly string[], mode: MafftScoringMode): string[] { + const rows = new Map(); let index = -1, sequence = ""; + const finish = () => { + if (index < 0) return; + if (rows.has(index)) throw new Error("MAFFT returned a duplicate sequence identifier."); + rows.set(index, sequence.toUpperCase()); + }; + for (const raw of source.split(/\r?\n/)) { + const line = raw.trim(); if (!line) continue; + if (line.startsWith(">")) { + finish(); index = Number(line.slice(1).trim()); sequence = ""; + if (!Number.isSafeInteger(index) || index < 0 || index >= expected.length) + throw new Error("MAFFT returned an unknown sequence identifier."); + } else { + if (index < 0) throw new Error("MAFFT returned sequence data before its first header."); + sequence += line.replace(/\s/g, ""); + } + } + finish(); + if (rows.size !== expected.length) throw new Error("MAFFT returned the wrong number of sequences."); + const aligned = expected.map((input, row) => { + const original = input.replaceAll("-", "").toUpperCase(), output = rows.get(row)!; + const normalized = mode === "amino-acid" ? original.replaceAll("*", "X") : original; + if (output.replaceAll("-", "") !== normalized) throw new Error(`MAFFT changed sequence ${row + 1}.`); + let cursor = 0; + return [...output].map((residue) => residue === "-" ? "-" : original[cursor++]).join(""); + }); + const width = aligned[0]?.length ?? 0; + if (!width || aligned.some((row) => row.length !== width)) throw new Error("MAFFT returned a non-rectangular alignment."); + return aligned; +} + +export function mafftProfileAddition(profile: readonly string[], aligned: readonly string[]): ProfileAddition { + if (aligned.length !== profile.length + 1) throw new Error("MAFFT profile addition returned the wrong number of rows."); + const sampleRows = aligned.slice(0, profile.length); + const withoutGapColumns = (rows: readonly string[]) => { + const columns = Array.from({ length: rows[0]?.length ?? 0 }, (_, column) => column) + .filter((column) => rows.some((row) => row[column] !== "-")); + return rows.map((row) => columns.map((column) => row[column]).join("")); + }; + const before = withoutGapColumns(profile), after = withoutGapColumns(sampleRows); + if (before.some((row, index) => row !== after[index])) + throw new Error("MAFFT reference addition changed the existing sample protein alignment."); + return { profileRows: sampleRows, sequence: aligned[profile.length] }; +} diff --git a/src/mafft-msa-runtime.ts b/src/mafft-msa-runtime.ts index 5ca5ee6..e5e90c9 100644 --- a/src/mafft-msa-runtime.ts +++ b/src/mafft-msa-runtime.ts @@ -1,3 +1,5 @@ +import { mafftProfileAddition } from "./mafft-msa-codec.ts"; +import type { MafftRequest } from "./mafft-msa-codec.ts"; export interface MafftProgress { detail: string } /** Run MAFFT 7.520 FFT-NS-2 in an isolated worker with growable WASM memory. */ @@ -6,8 +8,18 @@ export async function runMafftFftnsMsa( scoringMode: "literal" | "nucleotide" | "amino-acid" = "nucleotide", onProgress?: (progress: MafftProgress) => void, ): Promise { + if (signal?.aborted) throw new DOMException("MAFFT alignment cancelled.", "AbortError"); if (sequences.length < 2) return [...sequences]; - if (scoringMode === "amino-acid") throw new Error("The reference-panel MAFFT runner expects nucleotide sequences."); + return runMafft({ sequences: sequences.map(String), scoringMode }, signal, onProgress); +} + +/** Equivalent to mafft --amino --add reference.fa sample-alignment.fa. */ +export async function runMafftProfileAdd(profileRows: readonly string[], sequence: string, signal?: AbortSignal) { + const aligned = await runMafft({ sequences: [...profileRows, sequence], scoringMode: "amino-acid", profileRows: profileRows.length }, signal); + return mafftProfileAddition(profileRows, aligned); +} + +async function runMafft(request: MafftRequest, signal?: AbortSignal, onProgress?: (progress: MafftProgress) => void): Promise { const worker = new Worker(new URL("./mafft-msa-worker.ts", import.meta.url), { type: "module" }); return new Promise((resolve, reject) => { let settled = false; @@ -15,9 +27,9 @@ export async function runMafftFftnsMsa( if (settled) return false; settled = true; signal?.removeEventListener("abort", abort); worker.terminate(); return true; }; - const abort = () => { if (finish()) reject(new DOMException("MAFFT panel alignment cancelled.", "AbortError")); }; + const abort = () => { if (finish()) reject(new DOMException("MAFFT alignment cancelled.", "AbortError")); }; worker.onmessage = (event: MessageEvent<{ type: "result" | "progress" | "error"; result?: string[]; detail?: string; message?: string }>) => { - if (event.data.type === "progress") { onProgress?.({ detail: event.data.detail ?? "MAFFT is aligning candidate sequences" }); return; } + if (event.data.type === "progress") { onProgress?.({ detail: event.data.detail ?? "MAFFT is aligning sequences" }); return; } if (!finish()) return; if (event.data.type === "error") reject(new Error(event.data.message || "MAFFT failed.")); else resolve(event.data.result ?? []); @@ -25,6 +37,6 @@ export async function runMafftFftnsMsa( worker.onerror = (event) => { if (finish()) reject(new Error(event.message || "The MAFFT worker stopped unexpectedly.")); }; if (signal?.aborted) { abort(); return; } signal?.addEventListener("abort", abort, { once: true }); - worker.postMessage({ sequences: sequences.map(String) }); + worker.postMessage(request); }); } diff --git a/src/mafft-msa-worker.ts b/src/mafft-msa-worker.ts index 85fa190..4fe033b 100644 --- a/src/mafft-msa-worker.ts +++ b/src/mafft-msa-worker.ts @@ -1,5 +1,8 @@ /// +import { mafftArguments, numericMafftFasta, parseMafftAlignment } from "./mafft-msa-codec.ts"; +import type { MafftRequest } from "./mafft-msa-codec.ts"; + import mafftJavascriptUrl from "/biowasm/mafft/disttbfast.mjs?url"; import mafftWasmUrl from "/biowasm/mafft/disttbfast.wasm?url"; @@ -18,55 +21,12 @@ type MafftFactory = (options: { printErr(line: string): void; }) => Promise; -interface MafftRequest { sequences: string[] } - const encoder = new TextEncoder(); -function numericFasta(sequences: readonly string[]) { - return sequences.map((sequence, index) => `>${index}\n${sequence}\n`).join(""); -} - -function parseNumericAlignment(source: string, expected: readonly string[]) { - const rows = new Map(); - let index = -1, sequence = ""; - const finish = () => { - if (index < 0) return; - if (rows.has(index)) throw new Error("MAFFT returned a duplicate sequence identifier."); - rows.set(index, sequence.toUpperCase()); - }; - for (const raw of source.split(/\r?\n/)) { - const line = raw.trim(); - if (!line) continue; - if (line.startsWith(">")) { - finish(); index = Number(line.slice(1).trim()); sequence = ""; - if (!Number.isSafeInteger(index) || index < 0 || index >= expected.length) - throw new Error("MAFFT returned an unknown sequence identifier."); - } else { - if (index < 0) throw new Error("MAFFT returned sequence data before its first header."); - sequence += line.replace(/\s/g, ""); - } - } - finish(); - if (rows.size !== expected.length) throw new Error("MAFFT returned the wrong number of sequences."); - const aligned = expected.map((input, row) => { - const output = rows.get(row)!; - if (output.replaceAll("-", "") !== input.toUpperCase()) throw new Error(`MAFFT changed sequence ${row + 1}.`); - return output; - }); - const width = aligned[0]?.length ?? 0; - if (!width || aligned.some((row) => row.length !== width)) throw new Error("MAFFT returned a non-rectangular alignment."); - return aligned; -} - -const FFT_NS_2_ARGUMENTS = [ - "-q", "0", "-E", "2", "-V", "-1.53", "-s", "0.0", "-W", "6", "-O", - "-C", "0-0", "-D", "-b", "62", "-g", "0", "-f", "-1.53", "-Q", "100.0", - "-h", "0", "-F", "-X", "0.1", "-x", "-1", "-i", "/input.fa", -] as const; - self.addEventListener("message", (event: MessageEvent) => { void (async () => { const sequences = event.data.sequences.map((sequence) => sequence.toUpperCase()); + const { scoringMode, profileRows } = event.data; if (sequences.length < 2) { self.postMessage({ type: "result", result: sequences }); return; } if (sequences.some((sequence) => !sequence.length || /[^A-Z?*.-]/.test(sequence))) throw new Error("MAFFT input contains an unsupported symbol."); @@ -90,11 +50,11 @@ self.addEventListener("message", (event: MessageEvent) => { } }, }); - runtime.FS.writeFile("/input.fa", encoder.encode(numericFasta(sequences))); + runtime.FS.writeFile("/input.fa", encoder.encode(numericMafftFasta(sequences, scoringMode))); try { - const status = runtime.callMain([...FFT_NS_2_ARGUMENTS]); + const status = runtime.callMain(mafftArguments(scoringMode, profileRows)); if (status) throw new Error(`MAFFT exited with status ${status}.`); } finally { try { runtime.FS.unlink("/input.fa"); } catch { /* best effort */ } } - self.postMessage({ type: "result", result: parseNumericAlignment(stdout.join("\n"), sequences) }); + self.postMessage({ type: "result", result: parseMafftAlignment(stdout.join("\n"), sequences, scoringMode) }); })().catch((cause) => self.postMessage({ type: "error", message: cause instanceof Error ? cause.message : String(cause) })); }); diff --git a/src/postprocess.ts b/src/postprocess.ts index 514976c..a2c4c36 100644 --- a/src/postprocess.ts +++ b/src/postprocess.ts @@ -1,7 +1,7 @@ import { runAlivibeMsa } from "./alivibe-msa-runtime.ts"; import { runIndependentPanelFilter } from "./independent-panel-filter-runtime.ts"; -import type { PanelFilterResult } from "./independent-panel-filter.ts"; -import { runMafftFftnsMsa } from "./mafft-msa-runtime.ts"; +import type { PanelFilterResult, ProfileAddition } from "./independent-panel-filter.ts"; +import { runMafftFftnsMsa, runMafftProfileAdd } from "./mafft-msa-runtime.ts"; import { inspectAlignment, translateAlignedNucleotides } from "./alignment-utils.ts"; import { collapseAlignment } from "./collapse.ts"; import { extractAndScorePanel } from "./panel-profile.ts"; @@ -28,7 +28,15 @@ export type PanelFilterRunner = (sequences: readonly string[], panelRows: readon workers?: number, onProgress?: (progress: { completed: number; total: number }) => void) => Promise; export interface PostprocessProgress { fraction: number; detail: string } +export interface FunctionalAlignmentRunner { + msa: MsaRunner; + addToProfile: (profileRows: readonly string[], sequence: string, signal?: AbortSignal) => Promise; +} +const defaultFunctionalAlignment: FunctionalAlignmentRunner = { msa: runMafftFftnsMsa, addToProfile: runMafftProfileAdd }; + export interface PostprocessOptions { + /** Passing-protein export MSA and reference addition; defaults to MAFFT. */ + functionalAlignment?: FunctionalAlignmentRunner; /** Keep false when collapse is run as its own cancellable pipeline stage. */ collapse?: boolean; onCollapseProgress?: (progress: PostprocessProgress) => void; @@ -198,7 +206,7 @@ interface FunctionalOutcome { aa?: string; /** Exact-position identity to the clipped aligned reference, rounded to two decimal places. */ referenceMatch?: number; - /** Codon-aware alignment, clipped to the first/last reference residue. */ + /** Codon-aware export alignment preserving the reference-trimmed sequence. */ alignedNt?: string; alignedAa?: string; } @@ -294,6 +302,7 @@ function preserveProteinSequences(inputs: readonly string[], aligned: readonly s */ async function functionalFilterBatch( referenceName: string, reference: string, sequences: string[], threshold: number, runMsa: MsaRunner, signal?: AbortSignal, + exportAlignment: FunctionalAlignmentRunner = defaultFunctionalAlignment, ): Promise { const outcomes: FunctionalOutcome[] = Array(sequences.length), coding: Array<{ index: number; sequence: string }> = []; sequences.forEach((raw, index) => { @@ -335,7 +344,30 @@ async function functionalFilterBatch( outcomes[entry.index] = { passed: !reasons.length, reasons, nt: trimmed, aa, alignedNt: queryRegion, alignedAa: translateAlignedNucleotides(queryRegion, 0), referenceMatch: Number(rawRatio.toFixed(2)) }; }); - return { outcomes, referenceName: resolvedReferenceName, referenceNt: referenceRegion, referenceAa: translateAlignedNucleotides(referenceRegion, 0) }; + // Screening and its match scores remain independent of export layout. Build + // a fresh MSA of only the passing, reference-trimmed proteins: rejected ORFs + // and their gap columns must not determine the functional-pass alignment. + const passing = outcomes.filter((outcome) => outcome.passed); + if (!passing.length) return { outcomes, referenceName: resolvedReferenceName, + referenceNt: referenceRegion, referenceAa: translateAlignedNucleotides(referenceRegion, 0) }; + const proteins = passing.map((outcome) => outcome.aa!); + const profile = proteins.length > 1 + ? await runScalableMsa(proteins, exportAlignment.msa, signal, 0, "amino-acid") : proteins; + if (profile.length !== proteins.length || profile.some((row, index) => + row.length !== profile[0].length || degap(row) !== proteins[index])) + throw new Error("The functional export alignment did not preserve its input proteins."); + const aligned = await exportAlignment.addToProfile(profile, translate(referenceCoding), signal); + if (aligned.profileRows.length !== passing.length || degap(aligned.sequence) !== translate(referenceCoding) + || aligned.profileRows.some((row, index) => row.length !== aligned.sequence.length || degap(row) !== proteins[index])) + throw new Error("The functional export reference addition did not preserve its input proteins."); + // Keep all residues of the already-trimmed passing sequences. Clipping this + // second alignment again could silently remove terminal sample residues. + passing.forEach((outcome, index) => { + outcome.alignedAa = aligned.profileRows[index]; + outcome.alignedNt = backtranslate(outcome.alignedAa, outcome.nt!); + }); + return { outcomes, referenceName: resolvedReferenceName, + referenceNt: backtranslate(aligned.sequence, referenceCoding), referenceAa: aligned.sequence }; } export async function postprocess( @@ -459,7 +491,7 @@ export async function postprocess( referenceAlignments: Object.assign({}, ...completed.map((output) => output.referenceAlignments)), collapseGroups: Object.assign({}, ...completed.map((output) => output.collapseGroups)), functionalFilterErrors: Object.assign({}, ...completed.map((output) => output.functionalFilterErrors)), collapseSeconds: 0 }; - return options.collapse === false ? combined : collapsePostprocess(combined, config, signal, options.onCollapseProgress, runMsa); + return options.collapse === false ? combined : collapsePostprocess(combined, config, signal, options.onCollapseProgress, runMsa, undefined, options.functionalAlignment); } /** Run family-count-preserving haplotype collapse as its own resumable stage. */ @@ -468,6 +500,7 @@ export async function collapsePostprocess( onProgress?: (progress: PostprocessProgress) => void, runMsa: MsaRunner = runAlivibeMsa, sampleNames?: ReadonlySet, + functionalAlignment: FunctionalAlignmentRunner = defaultFunctionalAlignment, ): Promise { const started = performance.now(), alignments = { ...output.alignments }, referenceAlignments = { ...output.referenceAlignments }; const collapseGroups = { ...output.collapseGroups }, summaries = output.summaries.map((summary) => ({ ...summary })); @@ -502,7 +535,7 @@ export async function collapsePostprocess( let batch: FunctionalBatchOutcome; try { batch = await functionalFilterBatch(sample.functionalReferenceSequence.name, sample.functionalReferenceSequence.sequence, - collapsedRows.map((row) => row.sequence), sample.functionalMatchOverride ?? config.parameters.functionalMatchThreshold, runMsa, signal); + collapsedRows.map((row) => row.sequence), sample.functionalMatchOverride ?? config.parameters.functionalMatchThreshold, runMsa, signal, functionalAlignment); } catch (cause) { if (signal?.aborted || (cause instanceof DOMException && cause.name === "AbortError")) throw cause; const message = cause instanceof Error ? cause.message : String(cause); diff --git a/tests/fixtures/README.md b/tests/fixtures/README.md new file mode 100644 index 0000000..98cfeb3 --- /dev/null +++ b/tests/fixtures/README.md @@ -0,0 +1,6 @@ +`functional-env.fasta` contains the reference and three ungapped, trimmed passing +coding sequences (variants 1, 9, and 45) from donor_1_REN in the public nanopore +demo: https://github.com/MurrellGroup/PORPIDpipeline/blob/nanopore/raw-reads/demo.fastq.gz. +It covers a full-length Env and shorter ORFs, including divergent N termini, for +MAFFT sequence preservation, profile preservation, and terminal-stop regression +checks. Headers are shortened to keep the fixture independent of display labels. diff --git a/tests/fixtures/functional-env.fasta b/tests/fixtures/functional-env.fasta new file mode 100644 index 0000000..f754b26 --- /dev/null +++ b/tests/fixtures/functional-env.fasta @@ -0,0 +1,131 @@ +>reference +ATGAGAGTGAAGGAGAAATATCAGCACTTGTGGAGATGGGGGTGGAGATGGGGCACCATGCTCCTTGGGATGTTGATGAT +CTGTAGTGCTACAGAAAAATTGTGGGTCACAGTCTATTATGGGGTACCTGTGTGGAAGGAAGCAACCACCACTCTATTTT +GTGCATCAGATGCTAAAGCATATGATACAGAGGTACATAATGTTTGGGCCACACATGCCTGTGTACCCACAGACCCCAAC +CCACAAGAAGTAGTATTGGTAAATGTGACAGAAAATTTTAACATGTGGAAAAATGACATGGTAGAACAGATGCATGAGGA +TATAATCAGTTTATGGGATCAAAGCCTAAAGCCATGTGTAAAATTAACCCCACTCTGTGTTAGTTTAAAGTGCACTGATT +TGAAGAATGATACTAATACCAATAGTAGTAGCGGGAGAATGATAATGGAGAAAGGAGAGATAAAAAACTGCTCTTTCAAT +ATCAGCACAAGCATAAGAGGTAAGGTGCAGAAAGAATATGCATTTTTTTATAAACTTGATATAATACCAATAGATAATGA +TACTACCAGCTATAAGTTGACAAGTTGTAACACCTCAGTCATTACACAGGCCTGTCCAAAGGTATCCTTTGAGCCAATTC +CCATACATTATTGTGCCCCGGCTGGTTTTGCGATTCTAAAATGTAATAATAAGACGTTCAATGGAACAGGACCATGTACA +AATGTCAGCACAGTACAATGTACACATGGAATTAGGCCAGTAGTATCAACTCAACTGCTGTTAAATGGCAGTCTAGCAGA +AGAAGAGGTAGTAATTAGATCTGTCAATTTCACGGACAATGCTAAAACCATAATAGTACAGCTGAACACATCTGTAGAAA +TTAATTGTACAAGACCCAACAACAATACAAGAAAAAGAATCCGTATCCAGAGAGGACCAGGGAGAGCATTTGTTACAATA +GGAAAAATAGGAAATATGAGACAAGCACATTGTAACATTAGTAGAGCAAAATGGAATAACACTTTAAAACAGATAGCTAG +CAAATTAAGAGAACAATTTGGAAATAATAAAACAATAATCTTTAAGCAATCCTCAGGAGGGGACCCAGAAATTGTAACGC +ACAGTTTTAATTGTGGAGGGGAATTTTTCTACTGTAATTCAACACAACTGTTTAATAGTACTTGGTTTAATAGTACTTGG +AGTACTGAAGGGTCAAATAACACTGAAGGAAGTGACACAATCACCCTCCCATGCAGAATAAAACAAATTATAAACATGTG +GCAGAAAGTAGGAAAAGCAATGTATGCCCCTCCCATCAGTGGACAAATTAGATGTTCATCAAATATTACAGGGCTGCTAT +TAACAAGAGATGGTGGTAATAGCAACAATGAGTCCGAGATCTTCAGACCTGGAGGAGGAGATATGAGGGACAATTGGAGA +AGTGAATTATATAAATATAAAGTAGTAAAAATTGAACCATTAGGAGTAGCACCCACCAAGGCAAAGAGAAGAGTGGTGCA +GAGAGAAAAAAGAGCAGTGGGAATAGGAGCTTTGTTCCTTGGGTTCTTGGGAGCAGCAGGAAGCACTATGGGCGCAGCCT +CAATGACGCTGACGGTACAGGCCAGACAATTATTGTCTGGTATAGTGCAGCAGCAGAACAATTTGCTGAGGGCTATTGAG +GCGCAACAGCATCTGTTGCAACTCACAGTCTGGGGCATCAAGCAGCTCCAGGCAAGAATCCTGGCTGTGGAAAGATACCT +AAAGGATCAACAGCTCCTGGGGATTTGGGGTTGCTCTGGAAAACTCATTTGCACCACTGCTGTGCCTTGGAATGCTAGTT +GGAGTAATAAATCTCTGGAACAGATTTGGAATCACACGACCTGGATGGAGTGGGACAGAGAAATTAACAATTACACAAGC +TTAATACACTCCTTAATTGAAGAATCGCAAAACCAGCAAGAAAAGAATGAACAAGAATTATTGGAATTAGATAAATGGGC +AAGTTTGTGGAATTGGTTTAACATAACAAATTGGCTGTGGTATATAAAATTATTCATAATGATAGTAGGAGGCTTGGTAG +GTTTAAGAATAGTTTTTGCTGTACTTTCTATAGTGAATAGAGTTAGGCAGGGATATTCACCATTATCGTTTCAGACCCAC +CTCCCAACCCCGAGGGGACCCGACAGGCCCGAAGGAATAGAAGAAGAAGGTGGAGAGAGAGACAGAGACAGATCCATTCG +ATTAGTGAACGGATCCTTGGCACTTATCTGGGACGATCTGCGGAGCCTGTGCCTCTTCAGCTACCACCGCTTGAGAGACT +TACTCTTGATTGTAACGAGGATTGTGGAACTTCTGGGACGCAGGGGGTGGGAAGCCCTCAAATATTGGTGGAATCTCCTA +CAGTATTGGAGTCAGGAACTAAAGAATAGTGCTGTTAGCTTGCTCAATGCCACAGCCATAGCAGTAGCTGAGGGGACAGA +TAGGGTTATAGAAGTAGTACAAGGAGCTTGTAGAGCTATTCGCCACATACCTAGAAGAATAAGACAGGGCTTGGAAAGGA +TTTTGCTATAA +>variant_1 +ATGAGAGTGATGGGGATCAGGAAGAATTATCAGGACTGGTGGAAAGGGGGCATCTTGCTCCTTGGGATATTGATGATTAG +TAGTGCTACAGAACAATTGTGGGTCACAGTCTATTATGGGGTACCTGTGTGGAAAGAGGCAACCACTACTCTATTTTGTG +CCTCAGATGCTAAAGGATATGAAACAGAGGCACATAATGTTTGGGCCACACATGCCTGTGTGCCCACAGACCCCAGCCCA +CAAGAAATAAAAATGGAAAATGTGACAGAACATTTTAACATGTGGAAAAATTCCATGGTAGAACAAATGCATGAGGATAT +AATCAGTTTATGGGATCAAAGTCTAAAGCCATGTGTAAAACTAACCCCACTCTGTGTTACTTTAAATTGCATTAATACTA +CTAGTAGTGGTGTGAACGGGACAATAGAGGATGGGGACATGAAAAACTGCTCTTTCAATATAACCACAGACATGAAAGAT +AAGGTGCAGAGAGAATATGCATTTTTCAACAAACTTGATATAGTTCCAATGGATAGTAGTGGTAGTACTAATGAAACTAA +CTATAGGTTGATAAGTTGTAACACCTCAGCCATTACACAGGCCTGTCCAAAGATATCCTTTGAGCCAATTCCAATACATT +ATTGTGCCCCAGCAGGTTTTGCAATTATAAAGTGTAATAATAAGACATTCAATGGAACAGGACCATGTACAAATGTCAGC +ACAGTACAATGCACACATGGAATTAAGCCAGTAGTATCAACTCAACTACTGTTAAATGGCAGTCTAGCAGAAGAAGATAT +AATAATTAGATCTGAAAATCTCACAAATAATGCTAAAACCATAATAGTGCAGCTGAAGGACCCTGTAGAAATTACTTGTA +TAAGACCCAACAACAATACAAGAAAGAGTATAACTATGGGACCAGGAAGAGCATTCTATACAACAGACATAATAGGAAAT +ATAAGGCAAGCACATTGTAACATTAGTAGAGCAAATTGGACTGACACTCTAAAAAAGATAGCTAGCAAATTAAGAGAAAA +ATTTAATAAAACAATAATCTTTGATCAGCCTTCAGGAGGAGACCTAGAAATTGTAATGCACAGCTTTAATTGTGGAGGGG +AATTTTTCTACTGTAATACAACACAGCTGTTCAATAGTACATGGGATTGTACTGGAACTGCATGTACATTGAATGGAACT +GTTCTTAATGGGATAAATAACACCAATATTACTCTCCCATGCAGGATAAAACAAATTATAAACATGTGGCAAAAAGTAGG +AAAAGCAATGTACAACCCTCCCATCAAAGGAGTAATTCAATGTACATCAAATATTACAGGGATAATATTAACAAGCGATA +ATGGTGGGATTGGAAATAATGAAAATGAAACCTTCAGACCTGCAGGAGGAGATATGAGGGACAATTGGAGAAGTGAATTA +TATAAATATAAAGTAGTAAGAATTGAGCCATTAGGAATAGCACCCACCAAGGCAAAAAGGAGAGTGGTGCAGAGAGAAAA +AAGAGCAGTAACATTAGGAGCTATGTTCATTGGGTTCTTGGGAGCAGCAGGAAGCACTATGGGCGCAGCGTCAGTGACGC +TGACGGTACAGGCCAGACAATTGCTGTCTGGTATAGTGCAACAGCAGAGCAATCTGCTGAAGGCTATTGAGGCGCAACAG +CATCTGTTGCAACTCACAGTCTGGGGCATCAAACAGCTCCAGGCAAGAGTCCTGGCTATGGAAAGATACCTAAGGGACCA +ACAGCTCCTGGGGATTTGGGGTTGCTCTGGAAAACTCATCTGCACCACTGCCGTGCCTTGGAATATTACTTGGAGTAATA +AAAATCTGACTACAATTTGGGAAAACATGACCTGGATGGAGTGGGAAAAAGAAATTGACAATTACACAGAGTTGATATAC +AGCTTACTTGAAAATTCGCAAATACAACAAGAAAAGAATGAACAAGAATTATTGGCATTGGACCAATGGGCAAATTTGTG +GAGTTGGTTTGACATAACAAAATGGCTGTGGTATATAAAAATATTCATGATGATAGTAGGAGGCTTGATAGGTTTAAAAA +TAGTTTTTGCTGTACTTTCTATAGTAAATAGAGTTAGGCAGGGATACTCACCATTATCGTTGCAGATCCGACCCCCAGCC +CCGAGGGGACCCGACAGGCCCGGAGGAATCGAAGAAGACGGTGGAGAGAGAGACAGAGGCGGATCCGAACGATTAGCGAG +TGGATTCTTGACACTTATCTGGGTCGACCTGAGGGACCTGTGCCTCTTCAGCTACCGCCGCTTGACAGACTTACTCTCGA +TTGTAACGAGGACTGTGGAACTTCTGGGACGCAGGGGGTGGGAAATCCTCAAATATTGGTGGAATCTCCTACAATATTGG +AGTCAGGAACTCAAGAATAGTGCTACTAGCTTGCTTAACGCTACAGCCATAGCAGTAGGTGAGGGAACAGATAGGATCAT +AGAGGTATTACAAAGATTTGGTAGGGGCATTCTCCACATACCTAGAAGAATAAGACAGGGCTTGGAAAGGGCCTTGCTGT +AA +>variant_2 +ATGAAAAATGTGACAGAACATTTTAACATGTGGAAAAATTCCATGGTAGAACAAATGCATGAGGATATAATCAGTTTATG +GGATCAAAGTCTAAAGCCATGTGTAAAACTAACCCCACTCTGTGTTACTTTAAATTGCATTAATACTACTAGTAGTGGTG +TGAACGGGACAATAGAGGATGGGGACATGAAAAACTGCTCTTTCAATATAACCACAGACATGAAAAATAAGGTGCAGAGA +GAATATGCATTTTTCAACAAACTTGATATAGTTCCAATGGATAGTAGTGGTAGTACTAATGAAACTAACTATAGGTTGAT +AAGTTGTAACACCTCAGCCATTACACAGGCCTGTCCAAAGATATCCTTTGAGCCAATTCCAATACATTATTGTGCCCCAG +CAGGTTTTGCAATTATAAAGTGTAATAATAAGACATTCAATGGAACAGGACCATGTACAAATGTCAGCACAGTACAATGC +ACACATAGAATTAAGCCAGTAGTATCAACTCAACTACTGTTAAATGGCAGTCTAGCAGAAGAAGATATAATAATTAGATC +TGAAAATCTCACAAATAATGCTAAAACCATAATAGTGCAGCTGAAGGACCCTGTAGAAATTACTTGTATAAGACCCAACA +ACAATACAAGAAAGAGTATAACTATGGGACCAGAAAGAGCATTCTATACAACAGACATAATAGGAAATATAAGGCAAGCA +CATTGTAACATTAGTAGAGCAAATTGGACTGACACTCTAAAAAAGATAGCTAGCAAATTAAGAGAAAAATTTAATAAAAC +AATAATCTTTGATCAGCCTTCAGGAGGAGACCTAGAAATTGTAATGCACAGCTTTAATTGTGGAGGGGAATTTTTCTACT +GTAATACAACACAGCTGTTCAATAGTACATGGGATTGTACTGGAACTGCATGTACATTGAATGAAACTGTTCTTAATGGG +ATAAATAACACCAATATTACTCTCCCATGCAGGATAAAACAAATTATAAACATGTGGCAAAAAGTAGAAAAAGCAATGTA +CAACCCTCCCCTCAAAGGAGTAATTCAATGTACATCAAATATTACAGGGATAATATTAACAAGCGATAATGGTGGGATTG +GAAATAATGAAAATGAAACCTTCAGACCTGCAGGAGGAGATATGAGGGACAATTGGAGAAGTAAATTATATAAATATAAA +GTAATAAGAATTGAGCCATTAGGAATAGCACCCACCAAGGCAAAAAGGAGAGTGGTGCAAAGAGAAAAAAGAGCAGTAAC +ATTAGGAGCTATGTTCATTGGGTTCTTGGGAGCAGCAGGAAGCACTATGGGCGCAGCGTCAGTGACGCTGACGGTACAGG +CCAGACAATTGCTGTCTGGTATAGTGCAACAGCAGAGCAATCTGCTGAAGGCTATTGAGGCGCAACAGCATCTGTTGCAA +CTCACAGTCTGGGGCATCAAACAGCTCCAGGCAAGAGTCCTGGCTATGGAAAGATACCTAAGGGACCAACAGCTCCTGGG +GATTTGGGGTTGCTCTGGAAAACTCATCTGCACCACTGCCGTGCCTTGGAATATTACTTGGAGTAATAAAAATCTGACTA +CAATTTGGGAAAACATGACCTGGATGGAGTGGGAAAAAGAAATTGACAATTACACAGAGTTGATATACAGCTTACTTGAA +AATTCGCAAATACAACAAGAAAAGAATGAACAAGAATTATTGGCATTGGACCAATGGGCAAATTTGTGGAGTTGGTTTGA +CATAACAAAATGGCTGTGGTATATAAAAATATTCATGATGATAGTAGGAGGCTTGATAGGTTTAAAAATAGTTTTTGCTG +TACTTTCTATAGTAAATAGAGTTAGGCAGGGATACTCACCATTATCGTTGCAGATCCGACCCCCAGCCCCGAGGGGACCC +GACAGGCCCGGAGGAATCGAAGAAGACGGTGGAGAGAGAGACAAAGGCGGATCCGAACGATTAGCGAGTGGATTCTTGAC +ACTTATCTGGGTCGACCTGAGGGACCTGTGCCTCTTCAGCTACCGCCGCTTGACAGACTTACTCTCGATTGTAACGAGGA +CTGTGGAACTTCTGGGACGCAGGGGGTGGGAAATCCTCAAATATTGGTGGAATCTCCTACAATATTGGAGTCAGGAACTC +AAGAATAGTGCTACTAGCTTGCTTAACGCTACAGCCATAGCAGTAGGTGAGGAAACAGATAGGATCATAGAGGTATTACA +AAGATTTGGTAGGGGCATTCTCCACATACCTAGAAGAATAAGACAGGGCTTGAAAAGGGCCTTGCTGTAA +>variant_3 +ATGATTAGTAGTGCTACAGAACAATTGTGGGTCACAGTCTATTATGGGGTACCTGTGTGGAAAGAGGCAACCACTACTCT +ATTTTGTGCCTCAGATGCTAAAGGATATGAAACAGAGGCACATAATGTTTGGGCCACACATGCCTGTGTGCCCACAGACC +CCAGCCCACAAAAAATAAAAATGGAAAATGTGACAGAACATTTTAACATGTGGAAAAATTCCATGGTAGAACAAATGCAT +GAGGATATAATCAGTTTATGGGATCAAAGTCTAAAGCCATGTGTAAAACTAACCCCACTCTGTGTTACTTTAAATTGCAT +TAATACTACTAGTAGTGGTGTGAACGGGACAATAGAGGATGGGGACATGAAAAACTGCTCTTTCAATATAACCACAGACA +TAAAAGATAAGGTGCAGAGAGAATATGCATTTTTCAACAAACTTGATATAGTTCCAATGGATAGTAGTAGTAGTACTAAT +GAAACTAACTATAGGTTGATAAGTTGTAACACCTCAGCCATTACACAGGCCTGTCCAAAGATATCCTTTGAGCCAATTCC +AATACATTATTGTGCCCCAGCAGGTTTTGCAATTATAAAGTGTAATAATAAGACATTCAATGGAACAGGACCATGTACAA +ATGTCAGCACAGTACAATGCACACATGGAATTAAGCCAGTAGTATCAACTCAACTACTGTTAAATGGCAGTCTAGCAGAA +GAAGATATAATAATTAGATCTGAAAATCTCACAAATAATGCTAAAACCATAATAGTGCAGCTGAAGGACCCTGTAGAAAT +TACTTGTATAAGACCCAACAACAATACAAGAAAGAGTATAACTATGGGACCAGGAAGAGCATTCTATACAACAGACATAA +TAGGAAATATAAGGCAAGCACATTGTAACATTAGTAGAGCAAATTGGACTGACACTCTAAAAAAGATAGCTAGCAAATTA +AGAGAAAAATTTAATAAAACAATAATCTTTGATCAGCCTTCAGGAGGAGACCTAGAAATTGTAATGCACAGCTTTAATTG +TGGAGGGGAATTTTTCTACTGTAATACAACACAGCTGTTCAATAGTACATGGGATTGTACTGGAACTGCATGTACATTGA +ATGGAACTGTTCTTAATGGAATAAATAACACCAATATTACTCTCCCATGCAGGATAAAACAAATTATAAACATGTGGCAA +AAAGTAGGAAAAGCAATGTACAACCCTCCCATCAAAGGAGTAATTCAATGTACATCAAATATTACAGGGATAATATTAAC +AAGCGATAATGGTGGGATTGGAAATAATGAAAATGAAACCTTCAGACCTGCAGGAGGAGATATGAGGGACAATTGGAAAA +GTGAATTATATAAATATAAAGTAGTAAGAATTGAGCCATTAGGAATAGCACCCACCAAGGCAAAAAGGAGAGTGGTGCAG +AGAGAAAAAAGAGCAGTAACATTAGGAGCTATGTTCATTGGGTTCTTGGAAGCAGCAGGAAGCACTATGGGCGCAGCGTC +AGTGACGCTGACGGTACAGGCCAGACAATTGCTGTCTGGTATAGTGCAACAGCAGAGCAATCTGCTGAAGGCTATTGAGG +CGCAACAGCATCTGTTGCAACTCACAGTCTGGGGCATCAAACAGCTCCAGGCAAGAGTCCTGGCTATGGAAAGATACCTA +AGGGACCAACAGCTCCTGGGGATTTGGGGTTGCTCTGGAAAACTCATCTGCACCACTGCCGTGCCTTGGAATATTACTTG +GAGTAATAAAAATCTGACTACAATTTGGGAAAACATGACCTGGATGGAGTGGGAAAAAGAAATTGACAATTACACAGAGT +TGATATACAGCTTACTTAAAAATTCGCAAATACAACAAGAAAAGAATGAACAAGAATTATTGGCATTGGACCAATGGGCA +AATTTGTGGAGTTGGTTTGACATAACAAAATGGCTGTGGTATATAAAAATATTCATGATGATAGTAGGAGGCTTGATAGG +TTTAAAAATAGTTTTTGCTGTACTTTCTATAGTAAATAGAGTTAGGCAGGGATACTCACCATTATCGTTGCAGATCCGAC +CCCCAGCCCCGAGGGGACCCGACAGGCCCGGAGGAATCGAAGAAGACGGTGGAGAGAGAGACAGAGGCGGATCCGAACGA +TTAGCGAGTGGATTCTTGACACTTATCTGGGTCGACCTGAGGGACCTGTGCCTCTTCAGCTACCGCCGCTTGACAGACTT +ACTCTCGATTGTAACGAGGACTGTGGAACTTCTGGGACGCAGGGGGTGGAAAATCCTCAAATATTGGTGGAATCTCCTAC +AATATTGGAGTCAGGAACTCAAGAATAGTGCTACTAGCTTGCTTAACGCTACAGCCATAGCAGTAGGTGAGGGAACAGAT +AGGATCATAGAGGTATTACAAAGATTTGGTAGGGGCATTCTCCACATACCTAGAAGAATAAGACAGGGCTTGGAAAGGGC +CTTGCTGTAA diff --git a/tests/ts-core.test.mjs b/tests/ts-core.test.mjs index 6f5354b..89e4a46 100644 --- a/tests/ts-core.test.mjs +++ b/tests/ts-core.test.mjs @@ -8,7 +8,7 @@ import { renderToStaticMarkup } from "react-dom/server"; import { build } from "rolldown"; import { gunzipSync } from "fflate"; -import { inspectAlignment, summarizeAlignmentChanges, translateAlignmentFasta, validateCorrectedAlignment } from "../src/alignment-utils.ts"; +import { inspectAlignment, summarizeAlignmentChanges, translateAlignedNucleotides, translateAlignmentFasta, validateCorrectedAlignment } from "../src/alignment-utils.ts"; import { referenceDisplayColumns } from "../src/alignment-regions.ts"; import { ALIVIBE_BRIDGE_VERSION, ALIVIBE_SOURCE_REVISION, assertAlivibeInitialLoad, assertAlivibeRoundTripTarget, @@ -22,6 +22,8 @@ import { buildExportArchive, SAMPLE_EXPORT_KINDS } from "../src/export-archive.t import { nameMatchingSlot, referenceMappingRecords, referenceSlots } from "../src/input-mapping.ts"; import { mapParsimonyMutations } from "../src/phylo-mutations.ts"; import { collapsePostprocess, functionalFilter, postprocess } from "../src/postprocess.ts"; +import { createDirectMafftRunner } from "../cli-src/direct-mafft.mjs"; +import { mafftProfileAddition, parseMafftAlignment } from "../src/mafft-msa-codec.ts"; import { addSequenceToProfile, filterQueriesAgainstPanel } from "../src/independent-panel-filter.ts"; import { extractAndScorePanel } from "../src/panel-profile.ts"; import { functionalFilterStats, porpidCallStats, sampleOverviewStats } from "../src/report-stats.ts"; @@ -35,6 +37,8 @@ import { staticTreeHighlighterSvg } from "../src/static-tree-highlighter.ts"; import { PROCESSING_HUM_LEVEL, ProcessingHum } from "../src/processing-hum.ts"; import { buildStoredConsensusThresholdReview, prepareConsensusThresholdReanalysis, restoreUntouchedThresholdStatuses } from "../src/consensus-threshold-reanalysis.ts"; +const testFunctionalAlignment = (msa) => ({ msa, addToProfile: async (profile, sequence) => addSequenceToProfile(profile, sequence) }); + function spoolRecord(sample, hash, umi = "AACCGGTT", name = "read", sequence = "ACGT") { const encoder = new TextEncoder(), umiBytes = encoder.encode(umi), nameBytes = encoder.encode(name), sequenceBytes = encoder.encode(sequence); const body = 20 + umiBytes.length + nameBytes.length + sequenceBytes.length * 2, output = new Uint8Array(body + 4); @@ -557,7 +561,7 @@ test("functional alignment is codon-aware and clipped to the reference endpoints const consensus = [{ ...base.consensuses[0], sequence: query }]; const runner = async (sequences, _signal, _iterations, mode) => mode === "amino-acid" ? ["MKK*--", "MKKAA*"] : [`${sequences[0]}------`, sequences[1]]; - const output = await postprocess(consensus, [], config, undefined, runner, 1); + const output = await postprocess(consensus, [], config, undefined, runner, 1, undefined, { functionalAlignment: testFunctionalAlignment(runner) }); const functional = inspectAlignment(output.alignments["sample_1/functional-nucleotide"], 1); const functionalReference = inspectAlignment(output.referenceAlignments["sample_1/functional-nucleotide"], 1); assert.equal(functional.columns, query.length); assert.equal(functionalReference.columns, query.length); @@ -593,7 +597,7 @@ test("functional amino-acid alignment recovers a structurally truncated MSA row ]; const runner = async (sequences, _signal, _iterations, mode) => mode === "amino-acid" ? [sequences[0], `${sequences[1].slice(0, -1)}-`] : [...sequences]; - const output = await postprocess(consensuses, [], config, undefined, runner, 1); + const output = await postprocess(consensuses, [], config, undefined, runner, 1, undefined, { functionalAlignment: testFunctionalAlignment(async (sequences) => [...sequences]) }); const aligned = inspectAlignment(output.alignments["sample_1/functional-nucleotide"], 1).records; assert.equal(aligned.length, 3); assert.deepEqual(aligned.slice(1).map((row) => row.sequence.replaceAll("-", "")), ["ATGAAATAA", "ATGGGGTAA"]); assert.deepEqual(output.functionalFilterErrors, {}); @@ -611,10 +615,11 @@ test("an irrecoverable functional-filter error is isolated to one sample", async "good/uncollapsed-nucleotide": ">good_a\nATGAAATAA\n>good_b\nATGGGGTAA\n", }, referenceAlignments: {}, collapseGroups: {}, functionalFilterErrors: {}, collapseSeconds: 0 }; let aminoCalls = 0; - const output = await collapsePostprocess(input, config, undefined, undefined, async (sequences, _signal, _iterations, mode) => { + const runner = async (sequences, _signal, _iterations, mode) => { if (mode === "amino-acid" && aminoCalls++ === 0) throw new Error("deliberate irrecoverable test failure"); return [...sequences]; - }); + }; + const output = await collapsePostprocess(input, config, undefined, undefined, runner, undefined, testFunctionalAlignment(runner)); assert.match(output.functionalFilterErrors.bad, /deliberate irrecoverable/); assert.equal(output.functionalFilterErrors.good, undefined); assert.equal(output.alignments["bad/functional-nucleotide"], undefined); assert.match(output.alignments["good/functional-nucleotide"], /^>ref\n/); assert.equal(output.summaries.find((row) => row.sample === "bad").functionalPassed, undefined); @@ -632,6 +637,52 @@ test("functional reference profile-add preserves the existing sample MSA", () => assert.equal(added.sequence.replaceAll("-", ""), "MKA*"); }); +test("MAFFT protein profile addition preserves demo residues, column relationships, and stop codons", async () => { + const sequences = (await readFile("tests/fixtures/functional-env.fasta", "utf8")).trim().split(/^>/m).slice(1) + .map((record) => record.split(/\r?\n/).slice(1).join("")); + const proteins = sequences.map((sequence) => translateAlignedNucleotides(sequence)); + const runner = createDirectMafftRunner(resolve("public/biowasm/mafft/disttbfast.mjs"), resolve("public/biowasm/mafft/disttbfast.wasm")); + const profile = await runner(proteins.slice(1), undefined, 0, "amino-acid"); + const added = await runner.addToProfile(profile, proteins[0]); + assert.deepEqual(added.profileRows.map((row) => row.replaceAll("-", "")), proteins.slice(1)); + assert.equal(added.sequence.replaceAll("-", ""), proteins[0]); + assert(added.profileRows.every((row) => row.includes("*")), "MAFFT must not drop stop codons"); + const columns = Array.from({ length: added.sequence.length }, (_, index) => index) + .filter((column) => added.profileRows.some((row) => row[column] !== "-")); + assert.deepEqual(added.profileRows.map((row) => columns.map((column) => row[column]).join("")), profile); + assert.throws(() => mafftProfileAddition(profile, [added.profileRows[1], added.profileRows[0], added.profileRows[2], added.sequence]), /changed the existing/); + assert.throws(() => parseMafftAlignment(">0\nMKA\n", ["MKA*"], "amino-acid"), /changed sequence/); + const aborted = new AbortController(); aborted.abort(); + await assert.rejects(runner.addToProfile(profile, proteins[0], aborted.signal), /cancelled/); +}); + +test("functional export realigns only passes and preserves terminal residues without changing screening scores", async () => { + const base = resultBundle(), reference = "ATGAAAGCTTAA", rejected = "ATG" + "GGT".repeat(9) + "TAA"; + const input = { records: [], summaries: [{ ...base.summaries[0] }], alignments: { + "sample_1/uncollapsed-nucleotide": `>one\n${reference.padEnd(rejected.length, "-")}\n>two\n${"ATGAAGGCTTAA".padEnd(rejected.length, "-")}\n>bad\n${rejected}\n`, + }, referenceAlignments: {}, collapseGroups: {}, functionalFilterErrors: {}, collapseSeconds: 0 }; + const config = { ...base.config, parameters: { ...base.config.parameters, deterministicSeed: 1n, functionalMatchThreshold: .5 }, + samples: [{ ...base.config.samples[0], panelSequences: [], functionalReferenceSequence: { name: "ref", sequence: reference } }] }; + const screen = async (sequences) => sequences.map((row) => row.padEnd(Math.max(...sequences.map((s) => s.length)), "-")); + const calls = []; + const output = await collapsePostprocess(input, config, undefined, undefined, screen, undefined, { + msa: async (sequences, _signal, _iterations, mode) => { + calls.push({ sequences: [...sequences], mode }); return sequences.map((row) => `${row}-`); + }, + addToProfile: async (profileRows, sequence) => ({ profileRows: [...profileRows], sequence: `-${sequence}` }), + }); + assert.deepEqual(output.functionalFilterErrors, {}); + assert.equal(output.summaries[0].functionalPassed, 2); + assert.deepEqual(calls, [{ sequences: ["MKA*", "MKA*"], mode: "amino-acid" }]); + const exported = inspectAlignment(output.alignments["sample_1/functional-nucleotide"], 1).records; + assert.equal(exported[0].name, "ref"); assert.equal(exported[0].sequence, `---${reference}`); + assert.deepEqual(exported.slice(1).map((row) => row.sequence.replaceAll("-", "")), + output.collapseGroups.sample_1.filter((group) => group.functionalPass).map((group) => group.trimmedNt)); + assert(exported.slice(1).every((row) => row.sequence.startsWith("ATG")), "reference terminal gaps must not trigger a second clipping"); + assert(output.collapseGroups.sample_1.filter((group) => group.functionalPass).every((group) => group.referenceMatch >= .5)); + assert(exported.every((row) => ![...row.sequence.matchAll(/-+/g)].some((match) => match[0].length % 3))); +}); + test("static tree highlights only variants strictly above ten percent", () => { const bundle = { alignments: { "sample/nucleotide": ">sample_v1_9\nAAAA\n>sample_v2_1\nAAAT\n" }, trees: { "sample/nucleotide": "(sample_v1_9:0.01,sample_v2_1:0.02);" },