From 459c21eb4f5979fd6358d05edd45a51e4c3bdda2 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Sun, 4 Oct 2026 23:18:28 -0400 Subject: [PATCH 1/2] Separate dense SV discovery from complete-input candidate support --- CHANGELOG.md | 10 ++ docs/dense-sv-support-spec.md | 58 +++++++ docs/sv-rna.md | 53 ++++++ isovar/__init__.py | 2 +- isovar/cli/isovar_sv_rna.py | 6 + isovar/sv_rna.py | 23 ++- isovar/sv_rna_dense.py | 303 ++++++++++++++++++++++++++++++++++ isovar/sv_rna_orf_export.py | 3 + tests/test_sv_rna_dense.py | 245 +++++++++++++++++++++++++++ 9 files changed, 701 insertions(+), 2 deletions(-) create mode 100644 docs/dense-sv-support-spec.md create mode 100644 isovar/sv_rna_dense.py create mode 100644 tests/test_sv_rna_dense.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 57a8c6cf..32e1c9b5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,16 @@ Behavior changes that can alter results or break callers, by release. Patch releases that only fix bugs or add fixtures are omitted; see the [commit history](https://github.com/openvax/isovar/commits/master) for those. +## 1.42.2 + +- Add opt-in `isovar sv-rna --dense-support` / `dense_support=True`: scan + original alignments on disk, select event-first sequence/placement bundles + under the discovery cap, then count direct junction and exact full-path/ORF + witnesses across the complete eligible input. Original records, missing + qualities, cell/UMI ambiguity and signal ancestry remain auditable. Discovery + and ORF-start assessment limits remain explicit; complete candidate support + does not imply exhaustive discovery. Defaults are unchanged (#437). + ## 1.39.11 - Fragments with conflicting ref/alt and third-allele mate observations no diff --git a/docs/dense-sv-support-spec.md b/docs/dense-sv-support-spec.md new file mode 100644 index 00000000..7e55a36f --- /dev/null +++ b/docs/dense-sv-support-spec.md @@ -0,0 +1,58 @@ +# Dense SV discovery and support (#437) + +Add an opt-in `dense_support` API / `--dense-support` CLI mode. Scan the entire +original alignment input into a temporary SQLite spool, retaining original SAM +records and input-scoped segment identities. Group supplementary/alternative +records before selection. Search the requested breakpoint/exon/explicit regions; +prioritize event-bearing groups over ordinary context, then select deterministic +sequence/placement representatives under `max_records`. Duplicate abundance and +coordinate order must not hide a late rare event. Keep omitted discovery groups, +failed segment builds, query/path/extension limits explicit. + +Run the existing reconstruction engine on the bounded discovery BAM. Stream the +spooled original groups again to find direct junction/clip observations for the +fixed candidate paths. Reuse the existing placement, exact-interval, read-filter, +cell/UMI and signal-lineage policies. Count full-path and ORF witnesses from all +available originals, including one-read cells and missing QUAL. Keep discovery +votes separate from exhaustive direct/full-interval support; unrelated background +must never count as event support. Retain original witness records and query +intervals, and preserve ambiguous placements and junction-base alternatives. + +The complete-input support pass cannot make omitted hypotheses assessed, establish +independent molecules, or demonstrate translation. Ordinary intron-competition +assessment remains scoped to discovery records and must say so. Scratch disk scales +with input size; retained support witnesses and the JSON output can scale with +actual candidate support. The default bounded API and schema stay compatible. + +Verify rare late events, order/cap-independent support for a fixed candidate, +original SAM/QUAL preservation, conflicting mates, alternative paths, filters, +explicit omissions, CLI/export behavior, and the complete lint/test/CI gates. + +## Verification evidence + +The regression suite includes the original ATP5MG--KMT2A short-read fixture, +whose pinned peptide stays `MAQFVRNLVEKTPALVNG`, and synthetic split-read cases +with alternate sequences, reversed partner orientation, missing QUAL, conflicting +mate labels, repeated records and a visible signal parent outside search regions. +Original-read interval matching and the normal ORF exporter are reused. + +A local synthetic scale check (macOS, Python 3.11.12, pysam 0.24.1) placed +13 two-piece event segments after repeated 90-base donor background records. +With `max_records=2` and assembly disabled, every run selected two discovery +records and counted all 13 event fragments. Timings measure dense reconstruction; +peak RSS covers the Python process, including input creation and imports. + +| Background records | Reconstruction seconds | Peak RSS (MiB) | +|---:|---:|---:| +| 1,000 | 0.035 | 120.7 | +| 100,000 | 2.066 | 155.3 | +| 1,000,000 | 37.365 | 160.7 | + +These are synthetic background-scale checks, not throughput estimates for noisy +long reads or witness-heavy events. They do not replace a rerun of the Sid pilot. +The audit runner in draft PR #433 can opt into this mode when it is ready. + +Testing also reproduced the default collector's SAM-payload identity collision, +filed as [#438](https://github.com/openvax/isovar/issues/438). Dense mode raises on +a conflicting original instead of selecting an arbitrary payload. Identical +repeat records still deduplicate, and interrupted/failed spools are cleaned up. diff --git a/docs/sv-rna.md b/docs/sv-rna.md index 34870a09..c78b986f 100644 --- a/docs/sv-rna.md +++ b/docs/sv-rna.md @@ -464,3 +464,56 @@ Sources: [SAM](https://samtools.github.io/hts-specs/SAMv1.pdf) and [fusion reconstruction assessment](https://pmc.ncbi.nlm.nih.gov/articles/PMC6802306/), [NCBI translation tables](https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi), [Kozak 1986](https://pubmed.ncbi.nlm.nih.gov/3943125/). + +## Dense inputs: discovery and complete-input support + +`--dense-support` (Python: `dense_support=True`) scans the **entire original +alignment file**, including records outside the nominated regions. It groups +actual records by segment on disk, gives event-bearing junctions and breakpoint +clips first priority, and reserves representative bundles for distinct +sequence/placement classes before spending the discovery budget on repeated +copies. A late event need not fall within the first `--max-records` records. +Supplementary and alternative records of a selected segment stay together; +an SA tag alone still supplies no observed partner ([SAM specification](https://samtools.github.io/hts-specs/SAMv1.pdf)). + +```bash +isovar sv-rna --input event.json --bam original-regional.bam --output dense.json \ + --dense-support --sv-scratch-dir /path/to/scratch --max-records 50000 +``` + +The temporary SQLite spool and discovery BAM require scratch space proportional +to the whole input, and are removed on success or failure. Without +`--sv-scratch-dir`, the system temporary directory is used. Only the bounded +discovery records and retained candidate witnesses are held in memory; +witness-heavy results and their JSON can still be large. Use an original +regional BAM to avoid unnecessarily scanning a whole-genome alignment file. + +After bounded reconstruction, a separate pass visits all original eligible +segments to count support for the fixed candidate paths: + +- `junctions[].direct_support` counts the complete input's observed junctions, + retaining alternative junction-base sequences. A shared junction does not + establish the complete candidate sequence. +- `paths[].full_interval_support` requires one exact, placement-compatible + observation spanning the **entire path and its junction**. The exploratory + ORFs' `full_interval_support` uses the same rule for each ORF interval. + Partial context reads do not become complete witnesses. +- Original SAM records, query intervals, missing QUAL, cell/UMI evidence and + producer-aware signal ancestry remain available. Visible conflicting mate + labels remain unresolved, including mates outside the discovery regions. +- `sequence_evidence.voting_support` and seed counts remain discovery counts. + ORF start/inclusion assessments and competing annotated splices also remain + scoped to discovery records; their scope is recorded explicitly. +- `support_acquisition` describes the complete-input pass and its exclusions + and segment-build notes. `discovery` reports selected and omitted regional + records. `discovery_record_limit`, query/path/extension limits and any + `support_segment_grouping_limit` remain explicit. ORF exports retain both + discovery and support-pass metadata. + +Complete support for a returned candidate does not make hypothesis discovery +exhaustive. A `no_candidate_paths` result with discovery limits remains +unassessed evidence, not absence of an event or expression. Counts remain +input-scoped segments, fragments and labels, not independent molecules. + +Support-pass completion refers to the supplied alignment file. It cannot recover +original reads omitted during upstream acquisition, filtering or processing. diff --git a/isovar/__init__.py b/isovar/__init__.py index 8f056fad..2b7ea425 100644 --- a/isovar/__init__.py +++ b/isovar/__init__.py @@ -10,7 +10,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -__version__ = "1.42.1" +__version__ = "1.42.2" from .allele_interpretations import reconcile_allele_interpretations diff --git a/isovar/cli/isovar_sv_rna.py b/isovar/cli/isovar_sv_rna.py index e0de3e1f..aac34453 100644 --- a/isovar/cli/isovar_sv_rna.py +++ b/isovar/cli/isovar_sv_rna.py @@ -44,6 +44,11 @@ def make_parser(prog="isovar sv-rna"): help="Base/small-indel alternatives also need this fraction of the best (default: %(default)s)") group.add_argument("--max-records", type=int, default=SV_MAX_RECORDS, help="Alignment records examined; reaching the limit is reported (default: %(default)s)") + group.add_argument("--dense-support", action="store_true", + help="Scan the full original input on disk, prioritize diverse event seeds, and count " + "candidate support across all eligible original records") + group.add_argument("--sv-scratch-dir", type=Path, + help="Existing scratch directory for --dense-support (space scales with the input)") group.add_argument("--max-queries", type=int, default=SV_MAX_QUERIES, help="Region/mate/supplementary lookups; reaching the limit is reported (default: %(default)s)") group.add_argument("--max-paths", type=int, default=SV_MAX_PATHS, @@ -95,6 +100,7 @@ def run(args=None, prog=None): min_alternative_fragments=options.min_alternative_fragments, min_alternative_fraction=options.min_alternative_fraction, min_local_variant_fraction=options.min_local_variant_fraction, max_records=options.max_records, + dense_support=options.dense_support, scratch_dir=options.sv_scratch_dir, max_queries=options.max_queries, max_paths=options.max_paths, assemble=options.assemble, breakpoint_window=options.breakpoint_window, max_breakpoint_shift=options.max_breakpoint_shift, max_extension_segments=options.max_extension_segments, diff --git a/isovar/sv_rna.py b/isovar/sv_rna.py index 1dac838c..3c3c4746 100644 --- a/isovar/sv_rna.py +++ b/isovar/sv_rna.py @@ -1324,7 +1324,8 @@ def reconstruct_sv_rna(bam, *, event_id, reference_name, donor, acceptor, region inclusion_min_splice_anchor_bases=SV_INCLUSION_MIN_SPLICE_ANCHOR_BASES, inclusion_min_base_quality=SV_INCLUSION_MIN_BASE_QUALITY, inclusion_min_mapping_quality=SV_INCLUSION_MIN_MAPPING_QUALITY, - inclusion_mapq_255_is_unique=SV_INCLUSION_MAPQ_255_IS_UNIQUE): + inclusion_mapq_255_is_unique=SV_INCLUSION_MAPQ_255_IS_UNIQUE, + dense_support=False, scratch_dir=None): """Reconstruct RNA paths for one nominated, oriented DNA adjacency. Paths start at unannotated RNA junctions (and, when soft-clipped bases are @@ -1402,6 +1403,14 @@ def reconstruct_sv_rna(bam, *, event_id, reference_name, donor, acceptor, region in those gates. False treats it as unavailable, per the SAM specification; excluded 255 records are then reported as ``mapq_255_excluded_from_inclusion``. + dense_support : bool + Scan the complete original input into a disk spool, select bounded, + event-first discovery representatives, then count candidate junction + and full-interval witnesses across every original eligible segment. + Discovery and intron-competition limits remain explicit. + scratch_dir : str or pathlib.Path, optional + Existing directory for the dense-mode temporary spool and discovery + BAM. Scratch space scales with the complete input, not max_records. Returns ------- @@ -1410,6 +1419,11 @@ def reconstruct_sv_rna(bam, *, event_id, reference_name, donor, acceptor, region frame and event-linkage evidence, original records and limitations. This is not the validated supplied-fusion result consumed by Vaxrank. """ + options = locals().copy() + if not isinstance(dense_support, bool): + raise ValueError("dense_support must be True or False") + if scratch_dir is not None and not dense_support: + raise ValueError("scratch_dir requires dense_support") if not all((event_id, reference_name, sample_id, source)) or not isinstance(event_provenance, dict) \ or not event_provenance: raise ValueError("SV reconstruction requires event, reference, sample, source and event provenance") @@ -1444,6 +1458,13 @@ def reconstruct_sv_rna(bam, *, event_id, reference_name, donor, acceptor, region if len({(r.annotation, r.transcript_id) for r in references}) != len(references): raise ValueError("Duplicate reference transcript identity") collector = ReadCollector() if read_collector is None else read_collector + if dense_support: + from .sv_rna_dense import reconstruct_dense_sv_rna + + options.pop("bam") + options.update(references=references, regions=tuple(regions), read_collector=collector, + peptide_lengths=lengths) + return reconstruct_dense_sv_rna(bam, options) models = [_Model(r) for r in references] annotated = set().union(*(m.junctions() for m in models)) reference_peptides = {m.protein[i:i + k] for m in models if m.protein diff --git a/isovar/sv_rna_dense.py b/isovar/sv_rna_dense.py new file mode 100644 index 00000000..727551ca --- /dev/null +++ b/isovar/sv_rna_dense.py @@ -0,0 +1,303 @@ +"""Disk-backed, bounded discovery followed by complete-input SV support. + +Only candidate witnesses are retained in memory during support reporting. The +SQLite spool keeps unrelated background and arbitrarily distant pieces on disk. +""" + +from collections import Counter, defaultdict +from collections.abc import Mapping +from contextlib import closing +from hashlib import sha256 +from itertools import groupby +from pathlib import Path +import sqlite3 +from tempfile import TemporaryDirectory + +import pysam + +from .read_identity import segment_identity +from .read_lineage import ReadLineage +from .read_metadata import record_evidence +from .sv_rna import ( + _Event, _Model, _Observations, _ineligible, _join_class, _record_gaps, _record_id, + _path_result, reconstruct_sv_rna, +) +from .sv_rna_orfs import _witnesses +from .sv_rna_relations import LINKAGE_STATUSES + + +class _DiskGroups(Mapping): + """Eligible original segment records, retrieved on demand for ancestry.""" + + def __init__(self, connection, header): + self.connection, self.header = connection, header + + def __getitem__(self, identity): + rows = self.connection.execute( + 'SELECT sam FROM records WHERE rg=? AND name=? AND segment=? AND excluded IS NULL ORDER BY rid', + identity) + records = [pysam.AlignedSegment.fromstring(sam, self.header) for sam, in rows] + if not records: + raise KeyError(identity) + return records + + def __iter__(self): + yield from self.connection.execute( + 'SELECT DISTINCT rg, name, segment FROM records WHERE excluded IS NULL ORDER BY rg, name, segment') + + def __len__(self): + return self.connection.execute( + 'SELECT count(*) FROM (SELECT DISTINCT rg, name, segment FROM records WHERE excluded IS NULL)').fetchone()[0] + + def templates(self): + rows = self.connection.execute( + 'SELECT rg, name, segment, sam FROM records WHERE excluded IS NULL ORDER BY rg, name, segment, rid') + for _, template in groupby(rows, key=lambda row: row[:2]): + yield [pysam.AlignedSegment.fromstring(row[3], self.header) for row in template] + + +def _spool(bam, db, collector): + db.execute('CREATE TABLE records (rid TEXT PRIMARY KEY, rg TEXT, name TEXT, segment INTEGER, sam TEXT, excluded TEXT)') + scanned, excluded = 0, Counter() + bam.reset() + for read in bam.fetch(until_eof=True): + scanned += 1 + reason = 'unnamed_record' if read.query_name is None else _ineligible(read, collector) + if reason: + excluded[reason] += 1 + identity = segment_identity(read) + sam, rid = read.to_string(), _record_id(read) + try: + db.execute('INSERT INTO records VALUES (?, ?, ?, ?, ?, ?)', (rid, *identity, sam, reason)) + except sqlite3.IntegrityError: + previous = db.execute('SELECT sam FROM records WHERE rid=?', (rid,)).fetchone()[0] + if previous != sam: + raise ValueError('Conflicting original SAM records share SV record identity %s' % rid) + if scanned % 10000 == 0: + db.commit() + db.execute('CREATE INDEX segments ON records (rg, name, segment, rid)') + db.commit() + return scanned, dict(sorted(excluded.items())) + + +def _search_regions(bam, options): + lengths = dict(zip(bam.references, bam.lengths)) + regions = list(options['regions']) + [(r.contig, a, b) for r in options['references'] for a, b in r.exons] + for breakpoint in (options['donor'], options['acceptor']): + if breakpoint.contig not in lengths: + raise ValueError('Breakpoint contig %r is absent from the alignments' % breakpoint.contig) + regions.append((breakpoint.contig, max(0, breakpoint.position - options['breakpoint_window']), + min(lengths[breakpoint.contig], breakpoint.position + options['breakpoint_window']))) + spans = defaultdict(list) + for contig, start, end in regions: + if (contig not in lengths or type(start) is not int or type(end) is not int + or not 0 <= start < end <= lengths[contig]): + raise ValueError('Invalid or unavailable SV search interval: %r' % ((contig, start, end),)) + spans[contig].append((start, end)) + return spans + + +def _priority(records, collector, event, annotated): + """Event joins/clips first, other novel joins second, ordinary context last.""" + if len(records) > 1: + # Actual supplementary/alternative records stay together. An SA hint + # alone neither raises priority nor fabricates its missing record. + store = _Observations(records, collector, 1, {}, '', '') + for identity in store.groups: + for o in store.build(identity): + if event.clips(o.positions, len(o.sequence)): + return 0 + if any(_join_class(event, annotated, o.positions[i], o.positions[j], j - i - 1, kind)[0] + == 'event' for i, j, kind in o.breaks): + return 0 + priority = 2 + for read in records: + for left, right, kind, inserted in _record_gaps(read): + classes = [_join_class(event, annotated, left, right, len(inserted), kind), + _join_class(event, annotated, (right[0], right[1], '-'), + (left[0], left[1], '-'), len(inserted), kind)] + if any(use == 'event' for use, _ in classes): + return 0 + if not any(use in ('annotated', 'wobble') for use, _ in classes) and kind != 'D': + priority = 1 + if collector.use_soft_clipped_bases and any(op == 4 for op, _ in read.cigartuples): + if any(read.reference_name == b.contig and abs(edge - b.position) <= 1 + for b in (event.donor, event.acceptor) for edge in (read.reference_start, read.reference_end)): + return 0 + return priority + + +def _select(groups, db, options, collector, spans, event, annotated): + db.execute('CREATE TABLE discovery (rg TEXT, name TEXT, segment INTEGER, priority INTEGER, signature TEXT, size INTEGER)') + for template in groups.templates(): + members = defaultdict(list) + for read in template: + members[segment_identity(read)].append(read) + for identity, records in members.items(): + if not any(a < end and start < b for read in records for a, b in read.get_blocks() + for start, end in spans.get(read.reference_name, ())): + continue + # Ignore read name/labels, but retain library, sequence, placements, + # orientation, mate, quality and available path structure. + signature = sha256(repr((identity[0], identity[2], sorted( + (r.reference_name, r.reference_start, r.flag, r.cigarstring, r.query_sequence, + r.qual, r.mapping_quality) for r in records))).encode()).hexdigest() + db.execute('INSERT INTO discovery VALUES (?, ?, ?, ?, ?, ?)', + (*identity, _priority(records, collector, event, annotated), signature, len(records))) + db.commit() + # Reserve a threshold-sized bundle per sequence/placement class before + # spending the remaining budget on abundant copies of any one class. + query = '''SELECT rg, name, segment, size FROM ( + SELECT *, row_number() OVER (PARTITION BY priority, signature ORDER BY rg, name, segment) AS copy + FROM discovery) ORDER BY priority, (copy - 1) / ?, signature, copy''' + selected, count = [], 0 + for rg, name, segment, size in db.execute(query, (options['min_alternative_fragments'],)): + if count + size <= options['max_records']: + selected.extend(groups[rg, name, segment]) + count += size + total = db.execute('SELECT coalesce(sum(size), 0) FROM discovery').fetchone()[0] + return selected, total + + +def _positions(path): + positions = [None] * len(path['sequence']) + for block in path['blocks']: + bases = range(block['reference_start'], block['reference_end']) + if block['strand'] == '-': + bases = reversed(bases) + for q, base in zip(range(block['query_start'], block['query_end']), bases): + positions[q] = (block['contig'], base, block['strand']) + return tuple(positions) + + +def _support(groups, result, collector, event, annotated, options): + """Stream complete templates, keeping only candidate junction witnesses.""" + header = groups.header.to_dict() + joins = {(tuple(j['left']), tuple(j['right'])) for path in result['paths'] for j in path['junctions'] + if j['left'] is not None and j['right'] is not None} + placements = {(b['contig'], b['reference_start'], b['reference_end']) + for path in result['paths'] for b in path['blocks']} + records, reasons = [], Counter() + for template in groups.templates() if result['paths'] else (): + # No per-base path building for unrelated genomic background. + if not any(read.reference_name == contig and a < high and low < b + for read in template for a, b in read.get_blocks() for contig, low, high in placements): + continue + store = _Observations(template, collector, options['max_extension_segments'], header, + options['sample_id'], options['source']) + observations = [o for identity in store.groups for o in store.build(identity)] + reasons.update(store.reasons) + hit = any((o.positions[i], o.positions[j]) in joins + or event.relation(o.positions[i], o.positions[j], j - i - 1)[0] == 'breakpoint_junction' + for o in observations for i, j, _ in o.breaks) + hit |= any(event.clips(o.positions, len(o.sequence)) for o in observations) + if hit: + records.extend(template) # Include visible mates when resolving labels. + store = _Observations(records, collector, options['max_extension_segments'], header, + options['sample_id'], options['source']) + for identity in store.groups: + store.build(identity) + # Parents can be outside the candidates' intervals; keep original producer + # and ancestry resolution over the complete eligible input, without loading it. + store.lineage = ReadLineage(groups, header) + adjacency = [entry for (left, right), entries in store.joins.items() for entry in entries + if event.relation(left, right, entry[2] - entry[1] - 1)[0] == 'breakpoint_junction'] + clips = defaultdict(list) + for o in store.observations.values(): + for i, j, assignment in event.clips(o.positions, len(o.sequence)): + clips['donor' if assignment[0] == 0 else 'acceptor'].append((o, i, j, 'clip')) + used = {} + + def support(witnesses): + identities = {store.observations[w['observation']].identity for w in witnesses} + return dict(**store.cell_umi.support(identities), read_lineage=store.lineage.support(identities), + missing_quality_reads=len({store.observations[w['observation']].identity for w in witnesses + if store.observations[w['observation']].missing_qualities}), + scope='complete_input_direct_junction_observations', witnesses=witnesses) + + for path in result['paths']: + positions = _positions(path) + reported, evidence = _path_result(path['sequence'], positions, set(), store, event, annotated, adjacency, clips, + lambda *args: {}, lambda *args: None) + path['junctions'] = reported['junctions'] + used.update(evidence) + path['sequence_evidence']['voting_support_scope'] = 'bounded_discovery' + path['full_interval_support'] = support(_witnesses(path['sequence'], positions, 0, len(positions), + path['junctions'], evidence)) + for candidate in path['exploratory_orfs']['candidates']: + candidate['start_evidence']['observation_scope'] = 'bounded_discovery' + start, end = candidate['query_interval'] + crossed = [path['junctions'][i] for i in candidate['crossed_junctions']] + candidate['full_interval_support'] = support(_witnesses(path['sequence'], positions, start, end, + crossed, evidence)) + labels = {tuple(row['identity']): row for row in result['cell_umi_evidence']['reads']} + labels.update((tuple(row['identity']), row) for row in store.cell_umi.evidence()['reads']) + result['cell_umi_evidence']['reads'] = [labels[key] for key in sorted(labels)] + lineage = {tuple(row['identity']): row for row in result['read_lineage']['reads']} + lineage.update((tuple(row['identity']), row) for row in store.lineage.evidence()) + result['read_lineage'] = dict(scope='complete_eligible_input', reads=[lineage[key] for key in sorted(lineage)]) + for key, o in sorted(used.items()): + result['observations'][key] = dict(identity=list(o.identity), records=list(o.records), + reverse_complement=o.reverse_complement, + original_query_interval=list(o.query_interval), + missing_qualities=o.missing_qualities, secondary=o.secondary) + for read in store.groups[o.identity]: + rid = _record_id(read) + if rid in o.records: + result['original_records'][rid] = read.to_string() + result['record_evidence'][rid] = record_evidence(read) + return dict(reasons) + + +def reconstruct_dense_sv_rna(bam, options): + """Discover bounded hypotheses and recount their witnesses across the input.""" + options = dict(options) + scratch = options.pop('scratch_dir') + options['dense_support'] = False + collector = options['read_collector'] + spans = _search_regions(bam, options) + models = [_Model(r) for r in options['references']] + annotated = set().union(*(m.junctions() for m in models)) + event = _Event(options['donor'], options['acceptor'], options['max_breakpoint_shift'], annotated, + options['annotated_junction_tolerance']) + + with TemporaryDirectory(prefix='isovar-sv-', dir=scratch) as directory: + with closing(sqlite3.connect(str(Path(directory) / 'records.sqlite'))) as db: + db.execute('PRAGMA temp_store=FILE') + db.execute('PRAGMA cache_size=-8192') + scanned, excluded = _spool(bam, db, collector) + groups = _DiskGroups(db, bam.header) + selected, available = _select(groups, db, options, collector, spans, event, annotated) + unsorted, discovery = str(Path(directory) / 'discovery.unsorted.bam'), str(Path(directory) / 'discovery.bam') + with pysam.AlignmentFile(unsorted, 'wb', header=bam.header) as output: + for read in selected: + output.write(read) + pysam.sort('--no-PG', '-o', discovery, unsorted) + pysam.index(discovery) + with pysam.AlignmentFile(discovery) as discovery_bam: + result = reconstruct_sv_rna(discovery_bam, **options) + support_notes = _support(groups, result, collector, event, annotated, options) + result['parameters']['dense_support'] = True + result['acquisition']['input_scope'] = 'bounded_discovery_representatives' + result['support_acquisition'] = dict(scope='all_original_input_records', records_complete=True, + complete=not support_notes.get('segment_grouping_limit'), + records_scanned=scanned, unique_records=db.execute('SELECT count(*) FROM records').fetchone()[0], + eligible_segments=len(groups), + excluded_records=excluded, read_path_notes=support_notes) + result['discovery'] = dict(policy='event_first_sequence_placement_bundles.v1', + eligible_regional_records=available, selected_records=len(selected), + omitted_records=available - len(selected), + complete=available == len(selected) and not result['acquisition']['limitations']) + if available > len(selected): + result['limitations'] = sorted(set(result['limitations']) | {'discovery_record_limit'}) + if support_notes.get('segment_grouping_limit'): + result['limitations'] = sorted(set(result['limitations']) | {'support_segment_grouping_limit'}) + result['inclusion_competing_splice_scope'] = 'bounded_discovery' + result['seeds_scope'] = result['competing_annotated_junctions_scope'] = 'bounded_discovery' + result['observation_counts_scope'] = 'bounded_discovery' + rank = {status: i for i, status in enumerate(LINKAGE_STATUSES)} + result['paths'].sort(key=lambda path: ( + rank.get(path['event_linkage']['status'], len(rank)), + -max((j['direct_support']['fragments'] for j in path['junctions']), default=0), + path['sequence'], path['path_id'])) + return result diff --git a/isovar/sv_rna_orf_export.py b/isovar/sv_rna_orf_export.py index 1bbdf460..6777e36e 100644 --- a/isovar/sv_rna_orf_export.py +++ b/isovar/sv_rna_orf_export.py @@ -174,6 +174,9 @@ def export_sv_rna_orfs(result): export = {key: deepcopy(result[key]) for key in ( "event_id", "reference_name", "sample_id", "source", "event_provenance", "donor", "acceptor", "reference_models", "parameters", "limitations")} + for key in ("discovery", "support_acquisition", "inclusion_competing_splice_scope"): + if key in result: + export[key] = deepcopy(result[key]) export.update(schema="isovar.sv_rna_orfs.v4", reconstruction_schema=result["schema"], candidates=candidates, interval_convention="zero_based_half_open", evidence_identity_policy="sha256_of_domain_and_canonical_json; sample/source scoped; event excluded", diff --git a/tests/test_sv_rna_dense.py b/tests/test_sv_rna_dense.py new file mode 100644 index 00000000..b3d57a21 --- /dev/null +++ b/tests/test_sv_rna_dense.py @@ -0,0 +1,245 @@ +"""Complete-input support remains separate from bounded hypothesis discovery.""" + +from dataclasses import asdict +import json + +import pysam +import pytest + +from isovar import export_sv_rna_orfs +from isovar.cli.isovar_sv_rna import run as cli_run +from isovar.read_collector import ReadCollector +from tests.test_sv_rna import Scenario, aligned, record, write_bam + + +def dense_reads(scenario, n=13): + reads = [record('background-%04d' % i, '1', 2000, '90M', scenario.donor_ref.sequence[100:190]) + for i in range(250)] + for i in range(n): + pieces = aligned('event-%03d' % i, scenario.sequence[200:300], scenario.positions[200:300], + qualities=i != n - 1) + for read in pieces: + read.set_tag('CB', 'cell-%d' % i) + read.set_tag('UB', 'umi-%d' % i) + reads.extend(pieces) + return reads + + +def candidate(result): + path, = result['paths'] + return path + + +def test_late_rare_event_and_complete_support_do_not_depend_on_discovery_cap_or_order(tmp_path): + scenario = Scenario() + reads = dense_reads(scenario) + bam = tmp_path / 'dense.bam' + write_bam(bam, reads) + assert scenario.run(bam, max_records=2, assemble=False)['status'] == 'no_candidate_paths' + results = [] + for cap, ordering in [(2, reads), (4, list(reversed(reads))), (10, reads)]: + write_bam(bam, ordering) + result = scenario.run(bam, max_records=cap, assemble=False, dense_support=True, scratch_dir=tmp_path) + path = candidate(result) + assert path['sequence'] == scenario.sequence[200:300] + assert path['full_interval_support']['reads'] == path['full_interval_support']['fragments'] == 13 + assert path['full_interval_support']['cells'] == path['full_interval_support']['umis'] == 13 + assert path['full_interval_support']['missing_quality_reads'] == 1 + assert path['junctions'][0]['direct_support']['fragments'] == 13 + assert len(path['full_interval_support']['witnesses']) == 13 + assert path['sequence_evidence']['voting_support']['fragments'] <= cap // 2 + assert path['sequence_evidence']['voting_support_scope'] == 'bounded_discovery' + assert 'discovery_record_limit' in result['limitations'] + assert result['discovery']['omitted_records'] > 200 + assert result['support_acquisition']['complete'] + assert result['support_acquisition']['records_scanned'] == len(reads) + assert all('background' not in sam for sam in result['original_records'].values()) + assert {sam for sam in result['original_records'].values()} == { + read.to_string() for read in reads if read.query_name.startswith('event-')} + assert sorted(row['cell_barcode'] for row in result['cell_umi_evidence']['reads']) == [ + 'cell-%d' % i for i in sorted(range(13), key=lambda i: 'cell-%d' % i)] + results.append(path['full_interval_support']) + assert results[0] == results[1] == results[2] + # All spool/BAM files are cleaned on success. + assert not list(tmp_path.glob('isovar-sv-*')) + + +@pytest.mark.parametrize('strand', ['+', '-']) +def test_original_query_witnesses_and_orf_export_use_all_records(tmp_path, strand): + scenario = Scenario(strand) + reads = [] + for i in range(9): + reads += aligned('long-%d' % i, scenario.sequence, scenario.positions, qualities=False) + result = scenario.run(write_bam(tmp_path / 'long.bam', reads), dense_support=True, + max_records=4, assemble=False) + path = candidate(result) + assert path['full_interval_support']['fragments'] == 9 + assert path['full_interval_support']['missing_quality_reads'] == 9 + assert path['exploratory_orfs']['candidates'] + for orf in path['exploratory_orfs']['candidates']: + support = orf['full_interval_support'] + assert support['fragments'] == support['missing_quality_reads'] == 9 + assert len(support['witnesses']) == 9 + start, end = orf['query_interval'] + assert all(w['original_query_interval'] == [start, end] for w in support['witnesses']) + for orf in export_sv_rna_orfs(result)['candidates']: + assert orf['rna_support']['fragments'] == 9 + assert orf['rna_support']['missing_quality_reads'] == 9 + + +def test_filters_and_conflicting_mate_labels_are_preserved(tmp_path): + scenario = Scenario() + reads = [] + for i in range(5): + pieces = aligned('template-%d' % i, scenario.sequence[200:300], scenario.positions[200:300], flag=65) + for read in pieces: + read.set_tag('CB', 'cell-A') + read.set_tag('UB', 'umi-%d' % i) + reads += pieces + mate = record('template-%d' % i, '2', 15000, '80M', 'A' * 80, flag=129) + mate.set_tag('CB', 'cell-B') + mate.set_tag('UB', 'umi-%d' % i) + reads.append(mate) + # Duplicate and low-MAPQ event observations remain excluded in both passes. + reads += aligned('duplicate', scenario.sequence[200:300], scenario.positions[200:300], flag=1024) + low = aligned('low', scenario.sequence[200:300], scenario.positions[200:300]) + for read in low: + read.mapping_quality = 0 + reads += low + result = scenario.run(write_bam(tmp_path / 'mates.bam', reads), dense_support=True, + max_records=4, assemble=False) + support = candidate(result)['full_interval_support'] + assert support['reads'] == support['fragments'] == 5 + assert support['cells'] == support['umis'] == 0 + assert support['label_statuses'] == {'conflicting_template_labels': 5} + assert result['support_acquisition']['excluded_records'] + assert all('duplicate' not in sam and not sam.startswith('low\t') for sam in result['original_records'].values()) + + +def test_distinct_sequences_are_not_combined_into_one_supported_hypothesis(tmp_path): + scenario = Scenario() + original = scenario.sequence[200:300] + changed = original[:20] + ('A' if original[20] != 'A' else 'C') + original[21:] + reads = [read for i in range(5) for read in aligned('a-%d' % i, original, scenario.positions[200:300])] + reads += [read for i in range(4) for read in aligned('b-%d' % i, changed, scenario.positions[200:300])] + result = scenario.run(write_bam(tmp_path / 'alternatives.bam', reads), dense_support=True, + max_records=8, assemble=False) + assert {path['sequence']: path['full_interval_support']['fragments'] for path in result['paths']} == { + original: 5, changed: 4} + assert all(path['junctions'][0]['direct_support']['fragments'] == 9 for path in result['paths']) + + +def test_no_candidates_and_failed_scans_never_claim_negative_evidence(tmp_path): + scenario = Scenario() + reads = [record('ordinary', '1', 2000, '90M', scenario.donor_ref.sequence[100:190])] + result = scenario.run(write_bam(tmp_path / 'ordinary.bam', reads), dense_support=True, max_records=2) + assert result['status'] == 'no_candidate_paths' + assert result['discovery']['complete'] and result['support_acquisition']['complete'] + with pytest.raises(ValueError, match='scratch_dir requires'): + scenario.run(tmp_path / 'ordinary.bam', scratch_dir=tmp_path) + with pytest.raises(ValueError, match='positive integers'): + scenario.run(tmp_path / 'ordinary.bam', dense_support=True, max_records=0) + with pytest.raises(OSError): + scenario.run(tmp_path / 'ordinary.bam', dense_support=True, scratch_dir=tmp_path / 'absent') + + +def test_cli_dense_support_and_scratch_directory(tmp_path): + scenario = Scenario() + bam = write_bam(tmp_path / 'cli.bam', dense_reads(scenario, 5)) + inputs = dict(event_id='D--A', reference_name='test', sample_id='sample', + donor=asdict(scenario.donor), acceptor=asdict(scenario.acceptor), regions=[], + references=[asdict(scenario.donor_ref), asdict(scenario.acceptor_ref)], + event_provenance=dict(caller='synthetic DNA call')) + event, output = tmp_path / 'event.json', tmp_path / 'output.json' + event.write_text(json.dumps(inputs)) + cli_run(['--input', str(event), '--bam', str(bam), '--output', str(output), '--dense-support', + '--sv-scratch-dir', str(tmp_path), '--max-records', '2', '--no-assembly']) + result = json.loads(output.read_text()) + assert candidate(result)['full_interval_support']['fragments'] == 5 + assert result['parameters']['dense_support'] is True + + +def test_duplicate_fetches_and_alternative_placements_do_not_create_support(tmp_path): + scenario = Scenario() + reads = dense_reads(scenario, 4) + original = scenario.run(write_bam(tmp_path / 'original.bam', reads), max_records=4, + assemble=False, dense_support=True) + # Repeated identical records are one original alignment, not new molecules. + duplicate = scenario.run(write_bam(tmp_path / 'duplicate.bam', reads + reads), max_records=4, + assemble=False, dense_support=True) + assert candidate(original)['full_interval_support'] == candidate(duplicate)['full_interval_support'] + assert duplicate['support_acquisition']['records_scanned'] == 2 * len(reads) + assert duplicate['support_acquisition']['unique_records'] == len(reads) + + alternative = record('event-000', '1', 12000, '100M', scenario.sequence[200:300], flag=256) + ambiguous = scenario.run(write_bam(tmp_path / 'ambiguous.bam', reads + [alternative]), max_records=8, + assemble=False, dense_support=True) + assert candidate(ambiguous)['full_interval_support']['fragments'] == 4 + + +def test_conflicting_records_and_interrupted_scan_leave_no_partial_result(tmp_path, monkeypatch): + from isovar import sv_rna_dense + scenario = Scenario() + reads = dense_reads(scenario, 3) + bam = write_bam(tmp_path / 'interrupted.bam', reads) + original = sv_rna_dense._spool + + def interrupted(*args): + original(*args) + raise OSError('injected scan interruption') + + monkeypatch.setattr(sv_rna_dense, '_spool', interrupted) + with pytest.raises(OSError, match='injected scan interruption'): + scenario.run(bam, dense_support=True, scratch_dir=tmp_path) + assert not list(tmp_path.glob('isovar-sv-*')) + monkeypatch.setattr(sv_rna_dense, '_spool', original) + conflict = pysam.AlignedSegment.fromstring(reads[-1].to_string(), reads[-1].header) + conflict.set_tag('CB', 'another-cell') + bam = write_bam(tmp_path / 'conflicting.bam', reads + [conflict]) + with pytest.raises(ValueError, match='Conflicting original SAM'): + scenario.run(bam, dense_support=True, scratch_dir=tmp_path) + assert not list(tmp_path.glob('isovar-sv-*')) + + +def test_original_short_read_fixture_keeps_pinned_translation_and_support(tmp_path): + from tests.test_sv_rna import corpus_bam, run_corpus + _, bam, inputs = corpus_bam(tmp_path, 'ATP5MG--KMT2A') + original = run_corpus(bam, inputs) + result = run_corpus(bam, inputs, dense_support=True, max_records=10) + path = candidate(result) + assert path['sequence'] == candidate(original)['sequence'] + assert [t['amino_acids'] for t in path['translations']] == ['MAQFVRNLVEKTPALVNG'] + assert path['junctions'][0]['direct_support']['fragments'] == 2 + assert path['full_interval_support']['fragments'] <= 2 + assert result['support_acquisition']['records_complete'] + assert export_sv_rna_orfs(result)['support_acquisition'] == result['support_acquisition'] + + +def test_iterable_parameters_and_input_scoped_signal_parents(tmp_path): + from tests.test_sv_rna import HEADER + scenario = Scenario() + header = HEADER.to_dict() + header['RG'][0].update(PL='ONT', PG='dorado', LB='library', SM='sample') + header['PG'] = [dict(ID='dorado', PN='dorado', CL='dorado basecaller model input')] + header = pysam.AlignmentHeader.from_dict(header) + reads = [] + for i in range(3): + pieces = aligned('child-%d' % i, scenario.sequence[200:300], scenario.positions[200:300]) + for read in pieces: + read.set_tag('pi', 'parent') + read.set_tag('dx', 0) + read.set_tag('sp', i * 100) + reads += pieces + # The visible parent has invalid ancestry, beyond every candidate interval. + parent = record('parent', '2', 18000, '50M', 'A' * 50) + parent.set_tag('dx', 1) + reads.append(parent) + reads = [pysam.AlignedSegment.fromstring(r.to_string(), header) for r in reads] + bam = write_bam(tmp_path / 'parent.bam', reads, header=header) + result = scenario.run(bam, dense_support=True, max_records=2, assemble=False, + peptide_lengths=(n for n in (8, 9))) + support = candidate(result)['full_interval_support'] + assert support['fragments'] == 3 + assert support['read_lineage']['signal_groups'] == 0 + assert support['read_lineage']['statuses'] == {'unresolved_parent': 3} + assert any(r['identity'][1] == 'parent' for r in result['read_lineage']['reads']) From 222954d2413ba537460bea3822c34a4905c95e82 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Sun, 4 Oct 2026 23:33:22 -0400 Subject: [PATCH 2/2] Resolve original context metadata and prioritize placed event paths --- CHANGELOG.md | 2 +- docs/dense-sv-support-spec.md | 2 +- docs/sv-rna.md | 6 ++-- isovar/__init__.py | 2 +- isovar/sv_rna.py | 1 + isovar/sv_rna_dense.py | 58 ++++++++++++++++++++++++++--------- tests/test_sv_rna_dense.py | 40 ++++++++++++++++++++++++ 7 files changed, 91 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 32e1c9b5..ba184c2c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,7 @@ Behavior changes that can alter results or break callers, by release. Patch releases that only fix bugs or add fixtures are omitted; see the [commit history](https://github.com/openvax/isovar/commits/master) for those. -## 1.42.2 +## 1.44.0 - Add opt-in `isovar sv-rna --dense-support` / `dense_support=True`: scan original alignments on disk, select event-first sequence/placement bundles diff --git a/docs/dense-sv-support-spec.md b/docs/dense-sv-support-spec.md index 7e55a36f..09b05400 100644 --- a/docs/dense-sv-support-spec.md +++ b/docs/dense-sv-support-spec.md @@ -4,7 +4,7 @@ Add an opt-in `dense_support` API / `--dense-support` CLI mode. Scan the entire original alignment input into a temporary SQLite spool, retaining original SAM records and input-scoped segment identities. Group supplementary/alternative records before selection. Search the requested breakpoint/exon/explicit regions; -prioritize event-bearing groups over ordinary context, then select deterministic +prioritize placed event joins, then breakpoint clips, over ordinary context, then select deterministic sequence/placement representatives under `max_records`. Duplicate abundance and coordinate order must not hide a late rare event. Keep omitted discovery groups, failed segment builds, query/path/extension limits explicit. diff --git a/docs/sv-rna.md b/docs/sv-rna.md index c78b986f..43509c7d 100644 --- a/docs/sv-rna.md +++ b/docs/sv-rna.md @@ -469,8 +469,8 @@ Sources: [SAM](https://samtools.github.io/hts-specs/SAMv1.pdf) and `--dense-support` (Python: `dense_support=True`) scans the **entire original alignment file**, including records outside the nominated regions. It groups -actual records by segment on disk, gives event-bearing junctions and breakpoint -clips first priority, and reserves representative bundles for distinct +actual records by segment on disk, gives placed event junctions first priority, followed by breakpoint +clips, and reserves representative bundles for distinct sequence/placement classes before spending the discovery budget on repeated copies. A late event need not fall within the first `--max-records` records. Supplementary and alternative records of a selected segment stay together; @@ -502,6 +502,8 @@ segments to count support for the fixed candidate paths: producer-aware signal ancestry remain available. Visible conflicting mate labels remain unresolved, including mates outside the discovery regions. - `sequence_evidence.voting_support` and seed counts remain discovery counts. + `voting_observations` identifies the retained voters; their labels and ancestry + consult the complete input, including distant visible mates and parents. ORF start/inclusion assessments and competing annotated splices also remain scoped to discovery records; their scope is recorded explicitly. - `support_acquisition` describes the complete-input pass and its exclusions diff --git a/isovar/__init__.py b/isovar/__init__.py index 2b7ea425..c4f9ee0a 100644 --- a/isovar/__init__.py +++ b/isovar/__init__.py @@ -10,7 +10,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -__version__ = "1.42.2" +__version__ = "1.44.0" from .allele_interpretations import reconcile_allele_interpretations diff --git a/isovar/sv_rna.py b/isovar/sv_rna.py index 3c3c4746..6c55e5b1 100644 --- a/isovar/sv_rna.py +++ b/isovar/sv_rna.py @@ -1269,6 +1269,7 @@ def span(observation): somatic_causation_proven=False), sequence_evidence=dict( voting_support=store.cell_umi.support(segments), + voting_observations=sorted(voters), missing_quality_reads=len({o.identity for o in voting if o.missing_qualities}), secondary_reads=len({o.identity for o in voting if o.secondary}), assembly_phase="hypothesis_not_proven_long_range_phase"), diff --git a/isovar/sv_rna_dense.py b/isovar/sv_rna_dense.py index 727551ca..57c14789 100644 --- a/isovar/sv_rna_dense.py +++ b/isovar/sv_rna_dense.py @@ -15,9 +15,11 @@ import pysam +from .cell_umi import CellUmiEvidence from .read_identity import segment_identity from .read_lineage import ReadLineage from .read_metadata import record_evidence +from .rna_evidence import rna_support from .sv_rna import ( _Event, _Model, _Observations, _ineligible, _join_class, _record_gaps, _record_id, _path_result, reconstruct_sv_rna, @@ -55,6 +57,15 @@ def templates(self): for _, template in groupby(rows, key=lambda row: row[:2]): yield [pysam.AlignedSegment.fromstring(row[3], self.header) for row in template] + def template(self, fragment): + """All eligible original mates/placements of one consulted template.""" + groups = defaultdict(list) + for segment, sam in self.connection.execute( + 'SELECT segment, sam FROM records WHERE rg=? AND name=? AND excluded IS NULL ORDER BY segment, rid', + fragment): + groups[(*fragment, segment)].append(pysam.AlignedSegment.fromstring(sam, self.header)) + return groups + def _spool(bam, db, collector): db.execute('CREATE TABLE records (rid TEXT PRIMARY KEY, rg TEXT, name TEXT, segment INTEGER, sam TEXT, excluded TEXT)') @@ -98,19 +109,24 @@ def _search_regions(bam, options): def _priority(records, collector, event, annotated): - """Event joins/clips first, other novel joins second, ordinary context last.""" - if len(records) > 1: + """Placed event joins, breakpoint clips, other novel joins, then context.""" + clip = False + possible_clip = collector.use_soft_clipped_bases and any( + any(op == 4 for op, _ in read.cigartuples) + and any(read.reference_name == b.contig and abs(edge - b.position) <= 1 + for b in (event.donor, event.acceptor) for edge in (read.reference_start, read.reference_end)) + for read in records) + if len(records) > 1 or possible_clip: # Actual supplementary/alternative records stay together. An SA hint # alone neither raises priority nor fabricates its missing record. store = _Observations(records, collector, 1, {}, '', '') for identity in store.groups: for o in store.build(identity): - if event.clips(o.positions, len(o.sequence)): - return 0 + clip |= bool(event.clips(o.positions, len(o.sequence))) if any(_join_class(event, annotated, o.positions[i], o.positions[j], j - i - 1, kind)[0] == 'event' for i, j, kind in o.breaks): return 0 - priority = 2 + priority = 3 for read in records: for left, right, kind, inserted in _record_gaps(read): classes = [_join_class(event, annotated, left, right, len(inserted), kind), @@ -119,12 +135,8 @@ def _priority(records, collector, event, annotated): if any(use == 'event' for use, _ in classes): return 0 if not any(use in ('annotated', 'wobble') for use, _ in classes) and kind != 'D': - priority = 1 - if collector.use_soft_clipped_bases and any(op == 4 for op, _ in read.cigartuples): - if any(read.reference_name == b.contig and abs(edge - b.position) <= 1 - for b in (event.donor, event.acceptor) for edge in (read.reference_start, read.reference_end)): - return 0 - return priority + priority = 2 + return 1 if clip else priority def _select(groups, db, options, collector, spans, event, annotated): @@ -141,7 +153,8 @@ def _select(groups, db, options, collector, spans, event, annotated): # orientation, mate, quality and available path structure. signature = sha256(repr((identity[0], identity[2], sorted( (r.reference_name, r.reference_start, r.flag, r.cigarstring, r.query_sequence, - r.qual, r.mapping_quality) for r in records))).encode()).hexdigest() + r.qual, r.mapping_quality, r.get_tag('SA') if r.has_tag('SA') else '') + for r in records))).encode()).hexdigest() db.execute('INSERT INTO discovery VALUES (?, ?, ?, ?, ?, ?)', (*identity, _priority(records, collector, event, annotated), signature, len(records))) db.commit() @@ -230,11 +243,26 @@ def support(witnesses): crossed = [path['junctions'][i] for i in candidate['crossed_junctions']] candidate['full_interval_support'] = support(_witnesses(path['sequence'], positions, start, end, crossed, evidence)) - labels = {tuple(row['identity']): row for row in result['cell_umi_evidence']['reads']} + # Discovery-only context also needs original visible mates/parents. Its + # voting counts remain explicitly scoped to the bounded discovery set. + consulted = {tuple(row['identity']) for row in result['cell_umi_evidence']['reads']} + labels = {} + for fragment in sorted({identity[:2] for identity in consulted}): + original = CellUmiEvidence(groups.template(fragment), header, options['sample_id'], options['source']) + for identity in original.groups: + if identity in consulted: + original.segment(identity) + labels.update((tuple(row['identity']), row) for row in original.evidence()['reads']) labels.update((tuple(row['identity']), row) for row in store.cell_umi.evidence()['reads']) result['cell_umi_evidence']['reads'] = [labels[key] for key in sorted(labels)] - lineage = {tuple(row['identity']): row for row in result['read_lineage']['reads']} - lineage.update((tuple(row['identity']), row) for row in store.lineage.evidence()) + for path in result['paths']: + voting = {tuple(result['observations'][key]['identity']) + for key in path['sequence_evidence']['voting_observations']} + path['sequence_evidence']['voting_support'] = rna_support(voting, labels.__getitem__) + path['sequence_evidence']['label_resolution_scope'] = 'complete_eligible_input' + for row in result['read_lineage']['reads']: + store.lineage.segment(tuple(row['identity'])) + lineage = {tuple(row['identity']): row for row in store.lineage.evidence()} result['read_lineage'] = dict(scope='complete_eligible_input', reads=[lineage[key] for key in sorted(lineage)]) for key, o in sorted(used.items()): result['observations'][key] = dict(identity=list(o.identity), records=list(o.records), diff --git a/tests/test_sv_rna_dense.py b/tests/test_sv_rna_dense.py index b3d57a21..d30ee211 100644 --- a/tests/test_sv_rna_dense.py +++ b/tests/test_sv_rna_dense.py @@ -203,6 +203,7 @@ def interrupted(*args): def test_original_short_read_fixture_keeps_pinned_translation_and_support(tmp_path): from tests.test_sv_rna import corpus_bam, run_corpus + from isovar.read_end_inference import Adapter, ReadEndProfile _, bam, inputs = corpus_bam(tmp_path, 'ATP5MG--KMT2A') original = run_corpus(bam, inputs) result = run_corpus(bam, inputs, dense_support=True, max_records=10) @@ -213,6 +214,10 @@ def test_original_short_read_fixture_keeps_pinned_translation_and_support(tmp_pa assert path['full_interval_support']['fragments'] <= 2 assert result['support_acquisition']['records_complete'] assert export_sv_rna_orfs(result)['support_acquisition'] == result['support_acquisition'] + profile = ReadEndProfile('TruSeq', '1', (Adapter('TruSeq', 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC'),)) + trimmed = run_corpus(bam, inputs, dense_support=True, max_records=10, read_collector=ReadCollector( + use_soft_clipped_bases=True, read_end_profile=profile, trim_adapters=True)) + assert [p['sequence'] for p in trimmed['paths']] == [path['sequence']] def test_iterable_parameters_and_input_scoped_signal_parents(tmp_path): @@ -243,3 +248,38 @@ def test_iterable_parameters_and_input_scoped_signal_parents(tmp_path): assert support['read_lineage']['signal_groups'] == 0 assert support['read_lineage']['statuses'] == {'unresolved_parent': 3} assert any(r['identity'][1] == 'parent' for r in result['read_lineage']['reads']) + + +def test_discovery_context_labels_consult_original_mates_outside_search_regions(tmp_path): + scenario = Scenario() + reads = aligned('event', scenario.sequence[200:300], scenario.positions[200:300]) + context = record('context', '1', 2090, '60M', scenario.donor_ref.sequence[190:250], flag=65) + mate = record('context', '2', 18000, '60M', 'A' * 60, flag=129) + for read, cell in [(context, 'cell-A'), (mate, 'cell-B')]: + read.set_tag('CB', cell) + read.set_tag('UB', 'umi') + reads.extend([context, mate]) + result = scenario.run(write_bam(tmp_path / 'context.bam', reads), dense_support=True, + max_records=3, min_alternative_fragments=1, min_overlap=20) + rows = [row for row in result['cell_umi_evidence']['reads'] if row['identity'][1] == 'context'] + assert len(rows) == 2 + assert {row['status'] for row in rows} == {'conflicting_template_labels'} + assert all(path['sequence_evidence']['voting_support']['cells'] == 0 for path in result['paths']) + + +def test_placed_rare_event_precedes_unlinked_breakpoint_clips(tmp_path): + scenario = Scenario() + reads = [] + for i in range(40): + pieces = aligned('unlinked-%03d' % i, scenario.sequence[200:300], scenario.positions[200:300]) + for read in pieces: + read.set_tag('SA', '1,17000,+,100M,60,0;') + reads += pieces + reads += aligned('rare-linked', scenario.sequence[200:300], scenario.positions[200:300]) + result = scenario.run(write_bam(tmp_path / 'clips.bam', reads), *scenario.exact, + dense_support=True, max_records=2, assemble=False, + read_collector=ReadCollector(use_soft_clipped_bases=True)) + path = candidate(result) + assert path['sequence'] == scenario.sequence[200:300] + assert path['event_linkage']['status'] == 'breakpoint_junction' + assert path['full_interval_support']['fragments'] == 1