diff --git a/ChangeLog.md b/ChangeLog.md index d6a7ada..075b48e 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -6,6 +6,26 @@ Addresses [issue #4](https://github.com/BackofenLab/CopraRNA/issues/4), based on targeting `issue-1-python3-bugfixes`. This release completes the transition from the original 2.1.3 scripts to the Python-only CopraRNA2 runtime. +- **Correct genomic accessibility boundaries ([issue #6](https://github.com/BackofenLab/CopraRNA/issues/6)).** + Extend native target inputs on both sides by the target accessibility base-pair + span (`tAccL`, 100 nt by default), honoring shared `accL` and the native zero/window + semantics. Wrap circular genomes and stop at physical linear ends; preserve + strand, spliced CDS and short-truncated-UTR selection conventions. Restrict + interactions with `--tRegion` and native indexing so all reported positions, + seeds and structures refer to the requested sequence. Trim target CSV sequences + and normalize `E_norm` before filtering/statistics; reuse context for web exports. + Keep context sidecars in prediction signatures and retained job artifacts, and + invalidate old predictions. This intentionally changes accessibility energies and + potentially rankings; independently captured fixtures for unchanged extraction, + statistics and reports remain untouched. Test explicit boundary contracts and + real IntaRNAsTar interactions, including the effect of genomic folding context. +- **Simplify missing report fields and update DAVID access.** Following review of + [PR #7](https://github.com/BackofenLab/CopraRNA/pull/7) and the convention in + [upstream issue #22](https://github.com/PatrickRWright/CopraRNA/issues/22), write + unavailable values as empty cells in published result and web tables instead of + the internal `NA` marker. This makes CSV imports unambiguous without changing + statistical missing-value handling. Authenticate DAVID requests with the shared + email account instead of the former personal address. - **Remove static organism catalogs.** `CopraRNA_available_organisms.txt` supplied display names and an input allow-list; `kegg2refseqnew.csv` supplied assembly membership and arbitrary organism codes. Neither is necessary: resolve current diff --git a/README.md b/README.md index a7839ec..2624761 100644 --- a/README.md +++ b/README.md @@ -237,6 +237,32 @@ removed after success or failure. CopraRNA owns the native input/output paths, CSV format and required columns. It requests two interactions unless `outNumber` is supplied. Native options are not duplicated as CopraRNA flags. +Genomic targets include extra sequence on both sides for accessibility calculation: +the maximum target base-pair span `tAccL` (100 nt with IntaRNAsTar defaults), or the +shared `accL` setting. With `tAccL=0`, the flank uses `tAccW` (150 nt by default, +also configurable through `accW`). If both are zero, the context covers the whole +available linear replicon, or one full turn on each side of a circular target; +this unlimited setting can be expensive for genome-wide prediction. Circular +context wraps around the origin; linear context stops at the physical chromosome +ends. The same rules apply to 5′ UTR, 3′ UTR and spliced CDS targets, on both strands. +The existing exclusion of truncated UTR targets shorter than 150 nt is unchanged. + +CopraRNA sets IntaRNA's `--tRegion` to allow interactions only within the requested +region, and controls `tIdxPos0`/`qIdxPos0` so report coordinates start at 1. Targets +with different allowed intervals run in separate native batches. Accessibility is +computed on the extended sequence, including for web exports. FASTAs, interaction +CSVs, perfect-helix filtering, comparative statistics, alignments and plots use only +the requested sequence and coordinates; `E_norm` uses the requested target length. +Accessibility energies and probabilities retain the extra context's contribution, +so predicted interactions and rankings can change from earlier runs. + +Each target FASTA has a `.fa.context.json` sidecar containing the extended sequence +and its 0-based, half-open requested interval. These sidecars participate in cache +signatures and are retained in `IntaRNA/` after cleanup; temporary native inputs are +removed after success or failure. Existing predictions without this context must +be recomputed. Regenerate predictions before calling +`webserver.export_prediction()` independently on an old job. + Before fitting energy distributions, CopraRNA discards gap-free, stacking-only perfect helices spanning at least 90% of either RNA's effective length. Effective length is the smaller of sequence length and a positive `tAccW`/`qAccW`; zero @@ -298,18 +324,19 @@ coprarna run --workdir /data/jobs/balanced --srnaseq /data/my_sRNAs.fa \ Result columns include FDR-adjusted and comparative p-values, per-organism locus annotations and interaction coordinates, additional homologs and retained report fields. Coordinates follow IntaRNA's 1-based conventions relative to the extracted -sequences. A missing comparative score (`NA`) is not evidence of a significant +sequences. Final result and web tables write unavailable values as empty cells +rather than `NA`; a missing comparative score is not evidence of a significant interaction. Ranking identifies candidates for biological follow-up, not validated regulation. Review the 16S tree and weight warnings when an outlier dominates. All scientific diagrams are vector PDF/SVG. Extended PDFs retain every page, and the matching SVG stacks those pages. `--websrv` adds `coprarna_websrv_table.csv`, interaction details and SVG thumbnails. `--enrich N` writes optional functional -reports under `Enrichment/`. DAVID requires a reachable compatible SOAP service -and the account configured by the adapter; offline tests validate request -contracts, not live service availability. The optional interactive HTML enrichment -report loads upstream JavaScript when opened. Core prediction does not need a -browser or DAVID. +reports under `Enrichment/`. DAVID requires a reachable compatible SOAP service; +offline tests validate +request contracts, not live service availability. The optional interactive HTML +enrichment report loads upstream JavaScript when opened. Core prediction does not +need a browser or DAVID. ## Independent batches and helpers diff --git a/coprarna/conservation.py b/coprarna/conservation.py index 5f60725..33fc37d 100644 --- a/coprarna/conservation.py +++ b/coprarna/conservation.py @@ -6,7 +6,7 @@ import numpy as np from Bio import SeqIO -from .runtime import command, options, srna_input +from .runtime import command, srna_input from .rcompat import read_table, write_table, save_rdata, format_matrix_column from .rcompat import RNamedList from .statistics import r_number diff --git a/coprarna/david.py b/coprarna/david.py index 688d00d..98e02bd 100644 --- a/coprarna/david.py +++ b/coprarna/david.py @@ -5,6 +5,7 @@ WSDL = "https://david-d.ncifcrf.gov/webservice/services/DAVIDWebService?wsdl" ENDPOINT = "https://david-d.ncifcrf.gov/webservice/services/DAVIDWebService.DAVIDWebServiceHttpSoap11Endpoint/" +EMAIL = "rna@informatik.uni-freiburg.de" def _report(path, count, copra, client_factory=None): @@ -32,7 +33,7 @@ def emit(value): client = client_factory(WSDL) client.wsdl.services[0].setlocation(ENDPOINT) emit(client) - emit(client.service.authenticate("patrickrw@gmx.net")) + emit(client.service.authenticate(EMAIL)) emit( client.service.addList( ",".join(background[:count]), "ENTREZ_GENE_ID", "make_up", 0 diff --git a/coprarna/genome.py b/coprarna/genome.py index 36dec4d..be69207 100644 --- a/coprarna/genome.py +++ b/coprarna/genome.py @@ -1,6 +1,7 @@ """GenBank helpers preserving the legacy coordinate and boundary conventions.""" from pathlib import Path +from dataclasses import dataclass import re from .io import genbank_records, accession @@ -76,13 +77,51 @@ def sub(start, end): return f">{tag}\n{value}\n" -def parse_regions(path: Path, up: int, down: int, region: str) -> str: - """Extract each CDS target region from streamed GenBank records. +@dataclass(frozen=True) +class TargetRegion: + """A strand-oriented target with a 0-based, half-open prediction interval. - Coordinates follow the original 1-based inclusive flank conventions, with - strand, circular wrap and compound CDS handling. Each feature uses its own - record's sequence; the input file is not modified.""" - output = [] + ``sequence`` includes genomic accessibility context; ``requested`` excludes + it. Instances are in-memory values and never modify their GenBank source. + """ + + tag: str + sequence: str + start: int + end: int + + @property + def requested(self) -> str: + """Return only the requested target sequence, without accessibility flanks.""" + return self.sequence[self.start : self.end] + + +def _genomic_slice(sequence, start, end, circular): + # Context may cross either end or more than one circular origin. Linear + # records stop at their physical ends; never fabricate missing nucleotides. + if circular: + length = len(sequence) + offset = start % length + return (sequence * ((offset + end - start + length - 1) // length))[ + offset : offset + end - start + ] + return sequence[max(0, start) : min(len(sequence), end)] + + +def extract_targets( + path: Path, up: int, down: int, region: str, context: int | None = 0 +): + """Yield requested CDS regions with additional genomic accessibility context. + + Preserve the established requested-region sequence, including spliced CDS + order and omission of truncated UTRs shorter than 150 nt. Add ``context`` bases + at both biological ends, wrapping circular records and clipping linear ones. + None requests a full replicon length on each side for an unlimited span. + TargetRegion coordinates are 0-based, half-open within the extended sequence. + Stream compressed/plain GenBank without modifying it; reject negative context. + """ + if context is not None and context < 0: + raise ValueError("Accessibility context must be nonnegative") for record in genbank_records(path): source = str(record.seq) circular = record.annotations.get("topology") == "circular" @@ -106,10 +145,20 @@ def flank(pos, upstream, downstream, cds=False): if region == "5utr": pos = int(parts[0].start) + 1 if strand == 1 else int(parts[-1].end) - output.append(flank(pos, up, down)) + requested = flank(pos, up, down) + left, right = ( + (pos - up - 1, pos + down - 1) + if strand == 1 + else (pos + up, pos - down) + ) elif region == "3utr": pos = int(parts[-1].end) + 1 if strand == 1 else int(parts[0].start) - output.append(flank(pos, up, down)) + requested = flank(pos, up, down) + left, right = ( + (pos - up - 1, pos + down - 1) + if strand == 1 + else (pos + up, pos - down) + ) elif region == "cds": protein_region = str(feature.extract(record.seq)) # Perl spliced_seq is returned in biological order, independent of @@ -144,8 +193,45 @@ def flank(pos, upstream, downstream, cds=False): ).split("\n", 1)[1] if start == 1: after = after[1:] - output.append(before + protein_region + after) - return "".join(output) + requested = before + protein_region + after + left, right = ( + (int(parts[0].start) - up, int(parts[-1].end) + down) + if strand == 1 + else (int(parts[-1].end) + up, int(parts[0].start) - down) + ) + else: + raise ValueError(f"Unknown target region: {region}") + if not requested: + continue + sequence = requested.split("\n", 1)[1].rstrip("\n") + if not circular: + left, right = (min(len(source), max(0, p)) for p in (left, right)) + amount = len(source) if context is None else context + if strand == 1: + before = _genomic_slice(source, left - amount, left, circular) + after = _genomic_slice(source, right, right + amount, circular) + else: + before = _reverse_complement( + _genomic_slice(source, left, left + amount, circular) + ) + after = _reverse_complement( + _genomic_slice(source, right - amount, right, circular) + ) + yield TargetRegion( + tag, before + sequence + after, len(before), len(before) + len(sequence) + ) + + +def parse_regions(path: Path, up: int, down: int, region: str) -> str: + """Return requested CDS FASTA regions without accessibility context. + + Preserve the established 1-based flank conventions, strand, circular wrap, + spliced CDS and short-truncated-UTR exclusion. Stream compressed/plain + GenBank records without writing files; invalid coordinates raise ValueError. + """ + return "".join( + f">{r.tag}\n{r.requested}\n" for r in extract_targets(path, up, down, region) + ) def check_gene_cds(path: Path, display_name: str | None = None) -> str: diff --git a/coprarna/homology.py b/coprarna/homology.py index a9f026e..acb7bf5 100644 --- a/coprarna/homology.py +++ b/coprarna/homology.py @@ -102,6 +102,32 @@ def run_domclust(directory): (directory / "cluster.tab").write_text(table) +def _finalize_result_table(path): + """Publish a result CSV with empty missing cells and no internal sort column. + + Rewrite ``path`` in place, converting only complete ``NA`` fields so text that + happens to contain those letters is unchanged. Return the removed sort-column + values for the optional webserver alignment map. Reject empty or malformed + input rather than silently publishing a shifted table. + """ + lines = path.read_text().splitlines() + fields = [line.split(",") for line in lines] + if not fields or any(len(row) < 2 for row in fields): + raise ValueError(f"Malformed result table: {path}") + internal_sorting = [row[-1] for row in fields] + published = [ + ["" if cell == "NA" else cell for cell in row[:-1]] for row in fields + ] + text = "\n".join(",".join(row) for row in published) + "\n" + text = ( + text.replace(",Additional.homologs,", ",Additional homologs,") + .replace(",Amount.sampled", ",Amount sampled") + .replace("p.value", "p-value") + ) + path.write_text(text) + return internal_sorting + + def run_homology( srna: Path, up: int, @@ -319,18 +345,11 @@ def keep(line): plot_evolutionary_heatmap(directory) for filename in ("CopraRNA_result.csv", "CopraRNA_result_all.csv"): path = directory / filename - lines = path.read_text().splitlines() + internal_sorting = _finalize_result_table(path) if filename == "CopraRNA_result.csv" and int(opt["websrv"]): (directory / "CopraRNA_result.map_evo_align").write_text( - "\n".join(line.split(",")[-1] for line in lines) + "\n" + "\n".join(internal_sorting) + "\n" ) - text = "\n".join(",".join(line.split(",")[:-1]) for line in lines) + "\n" - text = ( - text.replace(",Additional.homologs,", ",Additional homologs,") - .replace(",Amount.sampled", ",Amount sampled") - .replace("p.value", "p-value") - ) - path.write_text(text) if int(opt["websrv"]): from .webserver import export_prediction diff --git a/coprarna/intarna.py b/coprarna/intarna.py index c8fea84..dce78dd 100644 --- a/coprarna/intarna.py +++ b/coprarna/intarna.py @@ -1,10 +1,12 @@ """IntaRNAsTar orchestration, compressed interactions and artifact filtering.""" from concurrent.futures import ThreadPoolExecutor +from dataclasses import asdict from pathlib import Path import csv import json import hashlib +import math import re import tempfile from Bio import SeqIO @@ -92,13 +94,17 @@ def filter_perfect_helices(path, windows=None): } -def prediction_signature(directory, target, query, *, native_version=None): +def prediction_signature( + directory, target, query, *, native_version=None, context=None +): """Hash scientific inputs and settings to prevent stale prediction reuse. Include tool version, personality, filter policy and file contents. Exclude absolute paths so equivalent jobs and inline/file settings have the same key. - Native callers supply the detected release instead of trusting stale options.""" - digest = hashlib.sha256(b"CopraRNA2:perfect-helix-filter=.9:v1") + Native callers supply the detected release instead of trusting stale options. + Include the context sidecar (extended sequences and allowed intervals) when + supplied. Only read inputs; no cache is published by this function.""" + digest = hashlib.sha256(b"CopraRNA2:perfect-helix-filter=.9:genomic-context:v2") version = ( native_version if native_version is not None @@ -107,6 +113,8 @@ def prediction_signature(directory, target, query, *, native_version=None): digest.update(version.encode()) digest.update(parameters.PERSONALITY.encode()) files = [Path(directory) / target, Path(directory) / query] + if context is not None: + files.append(Path(directory) / context) params = parameters.arguments(directory) if params: files.append(Path(params[1])) @@ -118,6 +126,100 @@ def prediction_signature(directory, target, query, *, native_version=None): return digest.hexdigest() +def _predict_targets(targets, query, output, directory, columns, *, runner=None): + """Predict extended TargetRegions and atomically publish requested-region CSV. + + Group targets by their 0-based allowed interval because IntaRNA 3 applies one + tRegion to a whole FASTA. Set tIdxPos0 so the requested sequence starts at 1 + (native negative indices skip zero), and restrict tRegion to 1..requested + length. Thus every native coordinate/seed/structure/bpList is already rebased. + Trim seq1 and normalize E_norm to requested lengths before any filtering or + downstream processing. Preserve native accessibility energies/probabilities. + Temporary FASTAs/CSVs are job-local and removed on failure; prior output is + replaced only after all groups validate. Callers retain their version checks. + """ + runner = runner or command + directory = Path(directory) + output = directory / output + groups = {} + seen = set() + for target in targets: + if target.tag in seen: + raise ValueError(f"Duplicate target locus tag: {target.tag}") + seen.add(target.tag) + if not 0 <= target.start < target.end <= len(target.sequence): + raise ValueError(f"Invalid prediction interval for {target.tag}") + groups.setdefault((target.start, target.end), []).append(target) + settings = parameters.values(directory) + with tempfile.TemporaryDirectory(prefix=".intarna-context-", dir=directory) as tmp: + temporary = Path(tmp) + combined = temporary / output.name + with open_text(combined, "wt") as stream: + writer = csv.DictWriter( + stream, columns.split(","), delimiter=";", lineterminator="\n" + ) + writer.writeheader() + for (start, end), records in groups.items(): + fasta = temporary / "targets.fa" + fasta.write_text("".join(f">{r.tag}\n{r.sequence}\n" for r in records)) + native = temporary / "native.csv.gz" + native.unlink(missing_ok=True) + runner( + [ + parameters.PERSONALITY, + *parameters.arguments(directory), + *([] if "outNumber" in settings else ["--outNumber=2"]), + "--target", + str(fasta.resolve()), + "--query", + str(query), + f"--tIdxPos0={-start if start else 1}", + "--qIdxPos0=1", + f"--tRegion=1-{end - start}", + "--outCsvCols=" + columns, + "--outMode=C", + "--out", + str(native.resolve()), + ], + directory, + ) + if not native.is_file(): + raise RuntimeError( + "IntaRNAsTar did not produce its interaction CSV" + ) + by_tag = {r.tag: r for r in records} + with open_text(native) as source: + reader = csv.DictReader(source, delimiter=";") + if reader.fieldnames != columns.split(","): + raise ValueError("Unexpected IntaRNA CSV columns") + for row in reader: + target = by_tag[row["id1"]] + if ( + not 1 + <= int(row["start1"]) + <= int(row["end1"]) + <= end - start + ): + raise ValueError( + f"IntaRNA interaction outside requested region: {target.tag}" + ) + if len(row["seq1"]) != len(target.sequence): + raise ValueError( + f"IntaRNA sequence length mismatch: {target.tag}" + ) + row["seq1"] = row["seq1"][start:end] + if "E_norm" in row: + # IntaRNA defines E_norm = E / ln(len(seq1)*len(seq2)). + # Keep native CSV's six significant-digit precision. + row["E_norm"] = format( + float(row["E"]) + / math.log(len(row["seq1"]) * len(row["seq2"])), + ".6g", + ) + writer.writerow(row) + combined.replace(output) + + def prepare_predictions( srna: Path, up: int, down: int, region: str, groups: list[str], directory: Path ) -> str: @@ -125,17 +227,29 @@ def prepare_predictions( Groups contain comma-separated replicon paths, with the primary chromosome first. Each organism runs independently using the same parameter file. Only - matching cached predictions are reused; filtering precedes all statistical work.""" + matching cached predictions are reused; filtering precedes all statistical work. + Write requested-region FASTAs for reports and .fa.context.json sidecars with + extended sequences/intervals for native prediction, cache keys and web export. + Extraction/native failures propagate without publishing a completed cache.""" native_version = parameters.check_version(directory) opt = options(directory) + context = parameters.target_context(directory) + target_regions = {} for group in groups: main = accession(group.split(",")[0]) target = directory / f"{main}_upfromstartpos_{up}_down_{down}.fa" - with target.open("w") as stream: - for filename in group.split(","): - stream.write( - genome.parse_regions(directory / filename, up, down, region) - ) + records = [ + r + for filename in group.split(",") + for r in genome.extract_targets( + directory / filename, up, down, region, context + ) + ] + target_regions[main] = records + target.write_text("".join(f">{r.tag}\n{r.requested}\n" for r in records)) + Path(str(target) + ".context.json").write_text( + json.dumps([asdict(r) for r in records]) + "\n" + ) queries = {} for record in SeqIO.parse(srna, "fasta"): queries[record.id.removeprefix("ncRNA_")] = record.id + ".fa" @@ -149,7 +263,11 @@ def predict(group): target = f"{main}_upfromstartpos_{up}_down_{down}.fa" output = target + ".intarna.csv.gz" signature = prediction_signature( - directory, target, queries[main], native_version=native_version + directory, + target, + queries[main], + native_version=native_version, + context=target + ".context.json", ) report_path = directory / (target + ".filter.json") try: @@ -160,25 +278,13 @@ def predict(group): not (directory / output).exists() or cached.get("native_signature") != signature ): - command( - [ - parameters.PERSONALITY, - *parameters.arguments(directory), - *([] if "outNumber" in windows else ["--outNumber", "2"]), - "--target", - target, - "--query", - queries[main], - "--outCsvCols", - tables.INTARNA_HEADER.replace(";", ","), - "--outMode=C", - "--out", - output, - ], + _predict_targets( + target_regions[main], + queries[main], + output, directory, + tables.INTARNA_HEADER.replace(";", ","), ) - if not (directory / output).is_file(): - raise RuntimeError(f"IntaRNAsTar did not produce {output}") report = filter_perfect_helices(directory / output, windows) report.update( native_signature=signature, diff --git a/coprarna/io.py b/coprarna/io.py index 3941dae..ad7cbb7 100644 --- a/coprarna/io.py +++ b/coprarna/io.py @@ -77,7 +77,7 @@ def genome_paths(directory): "aligned_sRNA.fa", "resolved_sRNA.fa", "ncRNA_*.fa", - "*_upfromstartpos_*_down_*.fa*", + "*_upfromstartpos_*_down_*.fa*", # Includes genomic-context JSON sidecars. "*.filter.json", "16s_sequences.*", "compatible.*", diff --git a/coprarna/parameters.py b/coprarna/parameters.py index 5930ddc..c66f8bb 100644 --- a/coprarna/parameters.py +++ b/coprarna/parameters.py @@ -80,8 +80,12 @@ def arguments(directory): def values(directory): - """Accessibility windows for the post-prediction artificial-helix check.""" - result = {"tAccW": "150", "qAccW": "150"} + """Read native settings with IntaRNAsTar accessibility defaults and aliases. + + No files are written. Resolve shared accW/accL settings as IntaRNA does and + reject conflicting per-RNA values rather than computing different flanks. + """ + result = {} configured = options(directory).get("intarna parameter file", "") path = ( Path(configured) if configured else Path(directory) / "IntaRNA_parameters.txt" @@ -94,4 +98,27 @@ def values(directory): if "=" in line: key, value = line.split("=", 1) result[key.strip()] = value.strip() + for shared, default in (("accW", "150"), ("accL", "100")): + for prefix in ("t", "q"): + specific = prefix + shared[0].upper() + shared[1:] + if shared in result: + if specific in result and int(result[specific]) != int(result[shared]): + raise ValueError(f"IntaRNA {shared} and {specific} must agree") + result[specific] = result[shared] + else: + result.setdefault(specific, default) return result + + +def target_context(directory): + """Return target accessibility flank length, or None for an unlimited span. + + IntaRNAsTar defaults to tAccL=100. Native zero means use tAccW (default 150); + both zero require full-replicon context. Read settings without modifying them + and reject invalid/negative spans before genomic extraction or native writes. + """ + settings = values(directory) + span, window = int(settings["tAccL"]), int(settings["tAccW"]) + if span < 0 or window < 0: + raise ValueError("IntaRNA accessibility span and window must be nonnegative") + return span or window or None diff --git a/coprarna/pipeline.py b/coprarna/pipeline.py index 4a442fc..9327dc6 100644 --- a/coprarna/pipeline.py +++ b/coprarna/pipeline.py @@ -141,7 +141,7 @@ def move_generated(path, destination): ): path.unlink() for subdir, patterns in [ - ("IntaRNA", ["*.fa.intarna.csv.gz", "*.filter.json"]), + ("IntaRNA", ["*.fa.intarna.csv.gz", "*.filter.json", "*.fa.context.json"]), ("Phylogeny", ["compatible.*", "16s_sequences.*"]), ("FASTA", ["*.fa"]), ("Regions_plots", ["*regions*"] + (["thumbnail_*"] if config.websrv else [])), diff --git a/coprarna/webserver.py b/coprarna/webserver.py index 81cfb8f..a0b92de 100644 --- a/coprarna/webserver.py +++ b/coprarna/webserver.py @@ -73,8 +73,10 @@ def shorten(value, isleft): def prepare_webserver_output(result: Path, interactions: Path, directory: Path) -> str: """Write the frontend table from comparative rows and native-derived duplexes. - Keep rows with missing comparative probabilities as NA; the more selective - sTar personality can leave valid organism interactions without a combined score.""" + Emit missing comparative probabilities as empty fields; the more selective + sTar personality can leave valid organism interactions without a combined score. + Write ``coprarna_internal_table.csv`` and raise ``ValueError`` for an empty input. + """ lines = result.read_text().splitlines() ui = interactions.read_text().splitlines() if len(lines) <= 1: @@ -103,7 +105,7 @@ def prepare_webserver_output(result: Path, interactions: Path, directory: Path) interaction = prune_interaction(*block[1:5], target, query) def number(value, precision): - return "NA" if value in ("", "NA") else format(float(value), precision) + return "" if value in ("", "NA") else format(float(value), precision) fields = [ str(rank), @@ -126,7 +128,7 @@ def number(value, precision): block[13], "0", ] - output.append(",".join(fields) + "\n") + output.append(",".join("" if field == "NA" else field for field in fields) + "\n") (directory / "coprarna_internal_table.csv").write_text("".join(output)) return "" @@ -191,14 +193,36 @@ def export_prediction(directory, ooi, up, down): Forward the same parameter file and interaction count as the primary run. Apply the same artificial-helix filter before choosing one retained native interaction per target for the frontend. This prevents a removed best hit - from reappearing when the comparative result used a retained suboptimal hit.""" + from reappearing when the comparative result used a retained suboptimal hit. + Require the main run's context sidecar, reuse its extended genomic sequences + and restrict/rebase interactions to the requested region. A missing or stale + sidecar fails before native prediction rather than losing accessibility context. + """ from Bio import SeqIO + from .genome import TargetRegion + from .intarna import _predict_targets parameters.check_version(directory) path = directory / f"{ooi}_upfromstartpos_{up}_down_{down}.fa" records = list(SeqIO.parse(path, "fasta")) lines = (directory / "CopraRNA_result.csv").read_text().splitlines()[1:] tags = [line.split(",")[2].split("(")[0] for line in lines if line.split(",")[2]] + context_path = Path(str(path) + ".context.json") + if not context_path.is_file(): + raise ValueError("Missing genomic accessibility context; rerun prepare-intarna") + context = { + r["tag"]: TargetRegion(**r) for r in json.loads(context_path.read_text()) + } + selected = [] + for tag in dict.fromkeys(t.casefold() for t in tags): + for record in records: + if tag == record.id.casefold(): + target = context[record.id] + if target.requested != str(record.seq): + raise ValueError( + f"Stale genomic accessibility context for {record.id}" + ) + selected.append(target) (directory / "CopraRNA_top_targets.fa").write_text( "".join( f">{r.id}\n{r.seq}\n" @@ -208,20 +232,13 @@ def export_prediction(directory, ooi, up, down): ) ) columns = "id1,seq1,id2,seq2,bpList,E,ED1,ED2,seedStart1,seedEnd1,seedStart2,seedEnd2,hybridDP,start1,end1,start2,end2" - command( - [ - parameters.PERSONALITY, - *parameters.arguments(directory), - "--target", - "CopraRNA_top_targets.fa", - "--query", - f"ncRNA_{ooi}.fa", - "--outMode=C", - "--outCsvCols=" + columns, - *([] if "outNumber" in parameters.values(directory) else ["--outNumber=2"]), - "--out=IntaRNA_web.csv.gz", - ], + _predict_targets( + selected, + f"ncRNA_{ooi}.fa", + "IntaRNA_web.csv.gz", directory, + columns, + runner=command, ) from .rcompat import read_table diff --git a/tests/test_accessibility_context.py b/tests/test_accessibility_context.py new file mode 100644 index 0000000..c773ef6 --- /dev/null +++ b/tests/test_accessibility_context.py @@ -0,0 +1,271 @@ +"""Requested-region contracts independent of the unchanged scientific oracles.""" + +import csv +from dataclasses import asdict +from io import StringIO +import json +import math +import os +import random + +from Bio import SeqIO +from Bio.Seq import Seq +from Bio.SeqFeature import CompoundLocation, SeqFeature, SimpleLocation +from Bio.SeqRecord import SeqRecord +import pytest + +from coprarna import genome, intarna, parameters, pipeline, tables +from coprarna.io import ensure_input_is_distinct, open_text, read_text +from coprarna.runtime import command + + +def settings(directory, contents=""): + path = directory / "native.cfg" + path.write_text(contents) + (directory / "CopraRNA_option_file.txt").write_text( + f"intarna parameter file:{path}\ncore count:1\n" + ) + return path + + +@pytest.mark.parametrize( + "contents, expected", + [ + ("", 100), + ("tAccL=60\n", 60), + ("accL=70\n", 70), + ("tAccL=0\n", 150), + ("tAccL=0\ntAccW=80\n", 80), + ("accL=0\naccW=90\n", 90), + ("accL=0\naccW=0\n", None), + ("qAccL=30\n", 100), + ("tAccL=60 # comment\ntAccW=80\n", 60), + ], +) +def test_context_uses_native_span_semantics(tmp_path, contents, expected): + path = settings(tmp_path, contents) + assert parameters.target_context(tmp_path) == expected + assert path.read_text() == contents + + +@pytest.mark.parametrize("contents", ["tAccL=-1\n", "accL=60\ntAccL=50\n"]) +def test_invalid_or_conflicting_context_fails(tmp_path, contents): + settings(tmp_path, contents) + with pytest.raises(ValueError): + parameters.target_context(tmp_path) + + +def genbank(directory, *, strand=1, circular=False, start=300, end=500, split=False): + rng = random.Random(42) + sequence = "".join(rng.choices("ACGT", k=1000)) + record = SeqRecord(Seq(sequence), id="NC_000001.1", name="NC_000001") + record.annotations = { + "molecule_type": "DNA", + "topology": "circular" if circular else "linear", + } + location = SimpleLocation(start, end, strand=strand) + if split: + parts = [ + SimpleLocation(start, start + 30, strand=strand), + SimpleLocation(end - 40, end, strand=strand), + ] + location = CompoundLocation(parts if strand == 1 else parts[::-1]) + record.features = [ + SeqFeature(location, type="CDS", qualifiers={"locus_tag": ["gene"]}) + ] + path = directory / "NC_000001.gb.gz" + with open_text(path, "wt") as stream: + SeqIO.write(record, stream, "genbank") + return path, sequence, location + + +@pytest.mark.parametrize("strand", [1, -1]) +@pytest.mark.parametrize("circular", [False, True]) +@pytest.mark.parametrize("region", ["5utr", "3utr", "cds"]) +@pytest.mark.parametrize( + "start,end", [(300, 500), (20, 220), (800, 980), (0, 200), (800, 1000)] +) +def test_genomic_context_strand_boundaries_and_requested_region( + tmp_path, strand, circular, region, start, end +): + path, source, location = genbank( + tmp_path, strand=strand, circular=circular, start=start, end=end + ) + up, down, padding = 100, 150, 60 + # Independent coordinates in the oriented genomic sequence. The UTR modes + # contain exactly up+down bases, with the same boundary as the public helper. + oriented = source if strand == 1 else str(Seq(source).reverse_complement()) + gene_start = start if strand == 1 else len(source) - end + gene_end = end if strand == 1 else len(source) - start + anchor = gene_end if region == "3utr" else gene_start + left = anchor - up + right = (gene_end if region == "cds" else anchor) + down + + def extract(a, b): + if circular: + return "".join(oriented[i % len(oriented)] for i in range(a, b)) + return oriented[max(0, a) : min(len(oriented), b)] + + requested = extract(left, right) + records = list(genome.extract_targets(path, up, down, region, padding)) + if not circular and region != "cds" and len(requested) < 150: + assert records == [] + return + (actual,) = records + assert actual.requested == requested + assert actual.sequence == extract(left - padding, right + padding) + assert actual.start == (padding if circular else min(padding, max(0, left))) + fasta = list( + SeqIO.parse(StringIO(genome.parse_regions(path, up, down, region)), "fasta") + ) + assert str(fasta[0].seq) == requested + + +@pytest.mark.parametrize("strand", [1, -1]) +def test_spliced_cds_keeps_exons_and_adds_external_context(tmp_path, strand): + path, source, location = genbank(tmp_path, strand=strand, split=True) + (actual,) = genome.extract_targets(path, 10, 20, "cds", 40) + if strand == 1: + before, after = source[250:300], source[500:560] + else: + before = str(Seq(source[500:550]).reverse_complement()) + after = str(Seq(source[240:300]).reverse_complement()) + assert actual.sequence == before + str(location.extract(Seq(source))) + after + assert (actual.start, actual.end) == (40, 140) + + +@pytest.mark.parametrize("circular", [False, True]) +def test_unlimited_span_and_multiple_circular_wraps(tmp_path, circular): + path, source, _ = genbank(tmp_path, circular=circular) + (actual,) = genome.extract_targets(path, 100, 150, "5utr", None) + if circular: + assert actual.start == 1000 + assert actual.sequence == "".join(source[i % 1000] for i in range(-800, 1450)) + else: + assert actual.sequence == source + assert (actual.start, actual.end) == (200, 450) + + +def test_context_and_interval_invalidate_cache(tmp_path): + settings(tmp_path) + (tmp_path / "targets.fa").write_text(">gene\nCCCCC\n") + (tmp_path / "query.fa").write_text(">query\nGGGGG\n") + context = tmp_path / "targets.fa.context.json" + signatures = [] + for sequence, start, end in [ + ("AAAAACCCCCAAAAA", 5, 10), + ("UUUUUCCCCCAAAAA", 5, 10), + ("UUUUUCCCCCAAAAA", 4, 9), + ]: + context.write_text( + json.dumps([asdict(genome.TargetRegion("gene", sequence, start, end))]) + ) + signatures.append( + intarna.prediction_signature( + tmp_path, + "targets.fa", + "query.fa", + native_version="IntaRNA 3.4.1", + context=context, + ) + ) + assert len(set(signatures)) == 3 + + +def test_context_sidecar_is_protected_and_cleanup_preserves_nested_input(tmp_path): + name = "NC_000001_upfromstartpos_200_down_100.fa.context.json" + original = tmp_path / "IntaRNA" / name + original.parent.mkdir() + original.write_text("original input\n") + generated = tmp_path / name + generated.write_text("generated context\n") + with pytest.raises(ValueError, match="collides"): + ensure_input_is_distinct(generated, tmp_path) + (tmp_path / "CopraRNA_option_file.txt").write_text(f"sRNA FASTA:{original}\n") + pipeline._cleanup(tmp_path, pipeline.Config()) + assert original.read_text() == "original input\n" + assert generated.read_text() == "generated context\n" + + +def test_failed_native_group_preserves_prior_csv_and_removes_temporaries(tmp_path): + settings(tmp_path) + output = tmp_path / "result.csv.gz" + output.write_bytes(b"prior result") + + def fail(*args): + raise RuntimeError("native failure") + + with pytest.raises(RuntimeError, match="native failure"): + intarna._predict_targets( + [genome.TargetRegion("gene", "AAAAACCCCCAAAAA", 5, 10)], + "query.fa", + output, + tmp_path, + tables.INTARNA_HEADER.replace(";", ","), + runner=fail, + ) + assert output.read_bytes() == b"prior result" + assert not list(tmp_path.glob(".intarna-context-*")) + + +@pytest.mark.skipif( + os.environ.get("COPRARNA_PIPELINE_TESTS") != "1", reason="opt-in native integration" +) +def test_native_context_changes_accessibility_and_preserves_report_coordinates( + tmp_path, +): + settings( + tmp_path, + "tAccW=80\ntAccL=60\nseedBP=3\noutMaxE=100\noutOverlap=B\noutNumber=3\n" + "tIdxPos0=500\nqIdxPos0=200\ntRegion=500-501\n", + ) + parameters.check_version(tmp_path) + # The C-rich requested target can fold with the added G-rich genomic flank. + target = genome.TargetRegion("gene", "G" * 6 + "AAAA" + "C" * 12 + "A" * 16, 10, 22) + short = genome.TargetRegion("short", "AAAA" + "C" * 8 + "AAAA", 4, 12) + (tmp_path / "query.fa").write_text(">query\nGGGGGGGGGGGG\n") + columns = tables.INTARNA_HEADER.replace(";", ",") + ",bpList" + calls = [] + + def record(args, directory): + calls.append(args) + return command(args, directory) + + intarna._predict_targets( + [target, short], "query.fa", "padded.csv.gz", tmp_path, columns, runner=record + ) + assert len(calls) == 2 + assert "--tRegion=1-12" in calls[0] and "--tIdxPos0=-10" in calls[0] + assert "--tRegion=1-8" in calls[1] and "--tIdxPos0=-4" in calls[1] + rows = list( + csv.DictReader(StringIO(read_text(tmp_path / "padded.csv.gz")), delimiter=";") + ) + assert {r["id1"] for r in rows} == {"gene", "short"} + for row in rows: + size = 12 if row["id1"] == "gene" else 8 + assert row["seq1"] == "C" * size + assert 1 <= int(row["start1"]) <= int(row["end1"]) <= size + for key in ("seedStart1", "seedEnd1"): + assert all(1 <= int(v) <= size for v in row[key].split(":")) + assert row["subseqDB"].startswith(row["start1"] + "C") + assert row["hybridDB"].startswith(row["start1"] + "|") + assert float(row["E_norm"]) == pytest.approx( + float(row["E"]) / math.log(size * 12), rel=1e-5 + ) + assert all( + 1 <= int(pair.split(",")[0][1:]) <= size + for pair in row["bpList"].split(":") + ) + intarna._predict_targets( + [genome.TargetRegion("gene", "C" * 12, 0, 12)], + "query.fa", + "bare.csv.gz", + tmp_path, + columns, + ) + bare = list( + csv.DictReader(StringIO(read_text(tmp_path / "bare.csv.gz")), delimiter=";") + ) + assert all(float(r["ED1"]) == 0 for r in bare) + assert any(float(r["ED1"]) > 0 for r in rows if r["id1"] == "gene") + assert not list(tmp_path.glob(".intarna-context-*")) diff --git a/tests/test_adapters.py b/tests/test_adapters.py index 4439cc5..a2fc4aa 100644 --- a/tests/test_adapters.py +++ b/tests/test_adapters.py @@ -13,10 +13,7 @@ @pytest.mark.parametrize("copra", [True, False]) -def test_david_requests_and_stdout_match_captured_source(copra, tmp_path, monkeypatch): - pytest.importorskip("suds.client") - import suds.client - +def test_david_requests_and_stdout_match_captured_source(copra, tmp_path): calls = [] class Service: @@ -42,7 +39,6 @@ def __init__(self, url): def __str__(self): return "offline DAVID client" - monkeypatch.setattr(suds.client, "Client", Client) path = tmp_path / "prediction.csv" if copra: path.write_bytes( @@ -59,6 +55,10 @@ def __str__(self): expected = json.loads( (FIXTURES / "david" / ("copra.json" if copra else "intarna.json")).read_text() ) + # The retained capture used the former personal DAVID account. Keep the + # independently captured request fixture intact and apply the reviewed account + # change explicitly at the contract boundary. + expected["calls"][2][1] = [david.EMAIL] actual = (david.coprarna_report if copra else david.intarna_report)(path, 5, Client) assert actual == expected["stdout"] assert json.loads(json.dumps(calls)) == expected["calls"] diff --git a/tests/test_reference_genomes.py b/tests/test_reference_genomes.py index 7fdf40f..168029e 100644 --- a/tests/test_reference_genomes.py +++ b/tests/test_reference_genomes.py @@ -62,3 +62,13 @@ def test_reference_genome_regions(path): assert parsed assert all(r.id in tags and len(r) > 0 for r in parsed) assert all(set(str(r.seq).upper()) <= set("ACGTRYSWKMBDHVN") for r in parsed) + extended = list(genome.extract_targets(path, 200, 100, region, 100)) + assert [(r.tag, r.requested) for r in extended] == [ + (r.id, str(r.seq)) for r in parsed + ] + assert all( + 0 <= r.start <= 100 and 0 <= len(r.sequence) - r.end <= 100 + for r in extended + ) + if record.annotations.get("topology") == "circular": + assert all(r.start == len(r.sequence) - r.end == 100 for r in extended) diff --git a/tests/test_review.py b/tests/test_review.py index 85cd030..759a79d 100644 --- a/tests/test_review.py +++ b/tests/test_review.py @@ -10,7 +10,7 @@ import pytest -from coprarna import genome, intarna, ncbi, parameters, pipeline, phylogeny +from coprarna import genome, homology, intarna, ncbi, parameters, pipeline, phylogeny from coprarna.cli import main from coprarna.helpers import HELPERS from coprarna.io import open_text, read_text, write_text @@ -304,7 +304,20 @@ def test_prediction_cache_signature_tracks_contents_not_paths(tmp_path): assert intarna.prediction_signature(tmp_path, "target.fa", "query.fa") != first -def test_report_retains_missing_comparative_pvalues(tmp_path): +def test_final_result_tables_publish_missing_values_as_empty_fields(tmp_path): + path = tmp_path / "CopraRNA_result.csv" + path.write_text( + "fdr,p.value,NC_000001,Annotation,Additional.homologs,Amount.sampled,initial_sorting\n" + "NA,NA,,containsNA,NA,0,7\n" + ) + assert homology._finalize_result_table(path) == ["initial_sorting", "7"] + assert path.read_text() == ( + "fdr,p-value,NC_000001,Annotation,Additional homologs,Amount sampled\n" + ",,,containsNA,,0\n" + ) + + +def test_web_report_publishes_missing_comparative_pvalues_as_empty_fields(tmp_path): from coprarna.webserver import prepare_webserver_output import shutil @@ -315,12 +328,16 @@ def test_report_retains_missing_comparative_pvalues(tmp_path): lines = result.read_text().splitlines() fields = lines[1].split(",") fields[:2] = ["NA", "NA"] + fields[-2] = "NA" lines[1] = ",".join(fields) result.write_text("\n".join(lines) + "\n") prepare_webserver_output(result, interactions, tmp_path) - assert (tmp_path / "coprarna_internal_table.csv").read_text().splitlines()[1].split( - "," - )[1:3] == ["NA", "NA"] + published = (tmp_path / "coprarna_internal_table.csv").read_text().splitlines()[ + 1 + ].split(",") + assert published[1:3] == ["", ""] + assert published[10] == "" + assert "NA" not in published def test_all_extended_plot_pages_are_exported_as_vectors(tmp_path, monkeypatch): @@ -352,6 +369,18 @@ def test_web_export_cannot_reintroduce_a_filtered_perfect_helix(tmp_path, monkey (tmp_path / "NC_000001_upfromstartpos_200_down_100.fa").write_text( ">tag\n" + "C" * 32 + "\n" ) + (tmp_path / "NC_000001_upfromstartpos_200_down_100.fa.context.json").write_text( + json.dumps( + [ + dict( + tag="tag", + sequence="A" * 100 + "C" * 32 + "A" * 100, + start=100, + end=132, + ) + ] + ) + ) (tmp_path / "CopraRNA_option_file.txt").write_text("intarna parameter file:\n") (tmp_path / "CopraRNA_result.csv").write_text( "fdr,p-value,NC_000001,Annotation,Additional homologs,Amount sampled\n" @@ -360,16 +389,19 @@ def test_web_export_cannot_reintroduce_a_filtered_perfect_helix(tmp_path, monkey def native(arguments, directory): assert "--outNumber=2" in arguments + assert "--tRegion=1-32" in arguments + assert "--tIdxPos0=-100" in arguments fields = next( arg.split("=", 1)[1] for arg in arguments if arg.startswith("--outCsvCols=") ).split(",") - with open_text(directory / "IntaRNA_web.csv.gz", "wt") as stream: + with open_text(arguments[arguments.index("--out") + 1], "wt") as stream: writer = csv.DictWriter(stream, fields, delimiter=";", lineterminator="\n") writer.writeheader() for pairs, score in [(32, -100), (10, -12)]: record = dict.fromkeys(fields, "0") record.update( row(pairs, query=32, target=32), + seq1="A" * 100 + "C" * 32 + "A" * 100, id1="tag", id2="query", bpList="".join(f"({i},{33 - i})" for i in range(1, pairs + 1)), diff --git a/tests/test_workflow.py b/tests/test_workflow.py index 2b83cb0..6911e29 100644 --- a/tests/test_workflow.py +++ b/tests/test_workflow.py @@ -182,10 +182,43 @@ def record(arguments, directory, **kwargs): for line in (job / "genome_paths.txt").read_text().splitlines() ) assert not list(job.glob("*.gb*")) + # Native context is retained separately; every downstream FASTA and CSV + # uses the originally requested 200+100 nt region and 1-based positions. + import csv + import json + from Bio import SeqIO + from coprarna.io import open_text + + sidecars = list(job.glob("*.fa.context.json")) + assert len(sidecars) == 3 + for sidecar in sidecars: + targets = json.loads(sidecar.read_text()) + requested = { + r.id: str(r.seq) + for r in SeqIO.parse( + str(sidecar).removesuffix(".context.json"), "fasta" + ) + } + assert all( + (r["start"], r["end"], len(r["sequence"])) == (60, 360, 420) + for r in targets + ) + assert all(r["sequence"][60:360] == requested[r["tag"]] for r in targets) + for output in [*job.glob("*.fa.intarna.csv.gz"), job / "IntaRNA_web.csv.gz"]: + with open_text(output) as stream: + rows = list(csv.DictReader(stream, delimiter=";")) + assert rows + assert all( + len(r["seq1"]) == 300 and 1 <= int(r["start1"]) <= int(r["end1"]) <= 300 + for r in rows + ) for pdf in job.glob("*.pdf"): assert pdf.with_suffix(".svg").is_file() assert calls and all(contents == parameter_text for _, _, contents in calls) assert len(calls) == 8 # Three organisms and one web export in each job. + assert all( + "--tIdxPos0=-60" in args and "--tRegion=1-300" in args for args, _, _ in calls + ) assert all( path.exists() if path == supplied else not path.exists() for _, path, _ in calls ) @@ -242,6 +275,13 @@ def test_corrected_workflow_builds_clusters_preserves_inputs_and_exports_vectors assert not (job / "cluster.members").exists() assert not (job / "target_alignments").exists() assert len(list((job / "IntaRNA").glob("*.csv.gz"))) == 3 + sidecars = list((job / "IntaRNA").glob("*.fa.context.json")) + assert len(sidecars) == 3 + assert all( + r["start"] == 100 and r["end"] == 400 and len(r["sequence"]) == 500 + for path in sidecars + for r in json.loads(path.read_text()) + ) tree = job / "Phylogeny/compatible.treefile" assert tree.is_file() and tree.with_suffix(".treefile.svg").is_file() for pdf in job.rglob("*.pdf"):