Fix genomic target accessibility boundaries

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

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
 

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

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

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:

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