Skip to content

Commit fb69548

Browse files
committed
.
1 parent 91c292d commit fb69548

2 files changed

Lines changed: 88 additions & 39 deletions

File tree

recmpox/diagnostic_snp.py

Lines changed: 27 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -471,16 +471,26 @@ def get_runs_and_breakpoints(
471471
positions_allegiances: List[Tuple[int, str]],
472472
diagnostic_snp_positions: List[int],
473473
min_consecutive: int = 1,
474+
ignore_other: bool = True,
474475
) -> Tuple[List[Tuple[int, int, str, int]], List[Tuple[int, int, str, str]]]:
475476
"""
476-
Build runs of consecutive ia/ib along diagnostic SNPs; "other" (ambiguous) ends a run.
477-
Only call a breakpoint when *both* the run before and the run after have >= min_consecutive SNPs.
478-
If min_consecutive is 1, all runs are considered (no consecutive-SNP filtering). SNPs classified
479-
as "other" are ignored for run boundaries.
477+
Build runs of consecutive ia/ib along diagnostic SNPs.
478+
479+
When ignore_other=True (default): "other" positions (N, gap, ambiguous base) are
480+
transparent to the run builder — a tract continues through them and only ends when
481+
an actual opposing-clade SNP is encountered. end_pos and n_snps reflect only the
482+
clade-matching positions; "other" positions inside a tract are excluded from its
483+
endpoints and SNP count. This makes tracts robust to patchy coverage: a stretch of
484+
N's inside an otherwise consistent Ia region will not split the tract.
485+
486+
When ignore_other=False: any "other" position ends the current run (original
487+
behaviour — more conservative, breaks tracts at every ambiguous site).
488+
489+
A breakpoint is only called when *both* flanking runs have >= min_consecutive SNPs.
480490
481491
Returns:
482492
runs: list of (start_pos, end_pos, clade, n_snps) with clade in ("ia", "ib").
483-
breakpoints: list of (pos_after_break, start_pos_next_run, clade_before, clade_after).
493+
breakpoints: list of (end_pos_before, start_pos_after, clade_before, clade_after).
484494
"""
485495
snp_positions = set(diagnostic_snp_positions)
486496
# Keep only SNP positions, sorted by position
@@ -494,17 +504,24 @@ def get_runs_and_breakpoints(
494504
if a not in ("ia", "ib"):
495505
i += 1
496506
continue
507+
current_clade = a
497508
start_pos = pos
509+
last_clade_pos = pos # last position that actually matched current_clade
498510
n_snps = 1
499511
i += 1
500512
while i < len(ordered):
501513
next_pos, next_a = ordered[i]
502-
if next_a != a:
514+
if next_a == current_clade:
515+
last_clade_pos = next_pos
516+
n_snps += 1
517+
i += 1
518+
elif next_a == "other" and ignore_other:
519+
# Transparent: skip without updating last_clade_pos or n_snps
520+
i += 1
521+
else:
522+
# Opposing clade (or "other" when ignore_other=False) → end this run
503523
break
504-
n_snps += 1
505-
i += 1
506-
end_pos = ordered[i - 1][0]
507-
runs.append((start_pos, end_pos, a, n_snps))
524+
runs.append((start_pos, last_clade_pos, current_clade, n_snps))
508525

509526
breakpoints: List[Tuple[int, int, str, str]] = []
510527
for j in range(len(runs) - 1):

recmpox/recmpox.py

Lines changed: 61 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -509,6 +509,29 @@ def _snp_positions_histogram_bins(
509509
return labels, counts
510510

511511

512+
def _genome_ruler_html(genome_length: int, min_width: int) -> str:
513+
"""Return an HTML ruler div with kbp tick marks proportional to genome_length."""
514+
if not genome_length:
515+
return ""
516+
# Pick a step size that gives 10–20 ticks
517+
step_bp = 10000
518+
for s in [1000, 2000, 5000, 10000, 20000, 50000, 100000, 200000, 500000]:
519+
if genome_length / s <= 20:
520+
step_bp = s
521+
break
522+
ticks = []
523+
pos = 0
524+
while pos <= genome_length:
525+
pct = pos / genome_length * 100
526+
label = "0" if pos == 0 else f"{pos // 1000}k"
527+
ticks.append(f'<span class="ruler-tick" style="left:{pct:.2f}%">{label}</span>')
528+
pos += step_bp
529+
# Always include a tick at the genome end if not already there
530+
if genome_length % step_bp != 0:
531+
ticks.append(f'<span class="ruler-tick" style="left:100%">{genome_length // 1000}k</span>')
532+
return f'<div class="strip-ruler" style="min-width:{min_width}px">{"".join(ticks)}</div>'
533+
534+
512535
def _write_results_html(
513536
out_path: Path,
514537
results: List[Dict[str, Any]],
@@ -649,29 +672,29 @@ def _ref_box(label: str, spec: Optional[str]) -> str:
649672

650673
# Diagnostic sites per sample: strip (genome position, color Ia/Ib/other) + table for ALL consensus genomes
651674
genome_length = results[0]["length"] if results else 0
675+
_first_alle = next((r.get("allegiances", []) for r in results if r.get("allegiances")), [])
676+
strip_min_w = max(600, len(_first_alle) * 2)
652677
rec_sites_html = ""
653-
n_segments_for_strip = 0
654678
for ri, r in enumerate(results):
655679
sample_id = r.get("id", "")
656680
allegiances = r.get("allegiances", [])
657681
rec_call = r.get("recombinant_call", "")
658682
if not allegiances:
659683
continue
660684
sorted_alle = sorted(allegiances, key=lambda x: x[0])
661-
if not n_segments_for_strip:
662-
n_segments_for_strip = len(sorted_alle)
663-
strip_min_w = max(400, n_segments_for_strip * 2)
664685
strip_segments = ""
665686
for (pos, allegiance) in sorted_alle:
666687
cls = "ia" if allegiance == "ia" else ("ib" if allegiance == "ib" else "other")
667688
lbl = ref1_label if allegiance == "ia" else (ref2_label if allegiance == "ib" else "other")
668-
strip_segments += f'<span class="strip-segment {cls}" title="{pos} {html_escape(lbl)}"></span>'
689+
pct = pos / genome_length * 100 if genome_length else 0
690+
strip_segments += f'<span class="strip-segment {cls}" title="{pos} bp – {html_escape(lbl)}" style="left:{pct:.3f}%"></span>'
669691
section_cls = "rec-sites-section" + (" recombinant" if rec_call == "potential recombinant" else "")
692+
ruler_html = _genome_ruler_html(genome_length, strip_min_w)
670693
rec_sites_html += (
671694
f'<div class="{section_cls}" data-row="{ri}" data-recombinant="{html_escape(rec_call)}">'
672695
f'<div class="rec-sites-row">'
673696
f'<span class="rec-sites-sample-id" title="{html_escape(sample_id)}">{html_escape(sample_id)}</span>'
674-
f'<div class="strip-cell"><div class="strip-genome" style="min-width:{strip_min_w}px" role="img" aria-label="Diagnostic sites along genome">{strip_segments}</div></div>'
697+
f'<div class="strip-cell"><div class="strip-genome" style="min-width:{strip_min_w}px" role="img" aria-label="Diagnostic sites along genome">{strip_segments}</div>{ruler_html}</div>'
675698
f'</div>'
676699
f'<details class="rec-sites-details"><summary>Show diagnostic site table (by tract)</summary>'
677700
f'<p class="threshold-note">Tracts = consecutive diagnostic sites with same classification. One row per tract.</p>'
@@ -722,7 +745,8 @@ def _ref_box(label: str, spec: Optional[str]) -> str:
722745
breakpoints: List[Tuple[int, int, str, str]] = []
723746
if allegiances:
724747
runs, breakpoints = get_runs_and_breakpoints(
725-
allegiances, diagnostic_snp_positions, min_consecutive=min_consecutive
748+
allegiances, diagnostic_snp_positions, min_consecutive=min_consecutive,
749+
ignore_other=True,
726750
)
727751
# Merge consecutive runs of the same clade (gaps = "other" ambiguous sites; treat as one tract)
728752
merged: List[Tuple[int, int, str, int]] = []
@@ -740,20 +764,16 @@ def _ref_box(label: str, spec: Optional[str]) -> str:
740764
merged_tracts[-1] = (merged_tracts[-1][0], end_pos, clade, merged_tracts[-1][3] + n_snps)
741765
else:
742766
merged_tracts.append((start_pos, end_pos, clade, n_snps))
767+
bp_strip_min_w = max(600, genome_length // 150) if genome_length else 600
743768
strip_segments = ""
744769
for j, (start_pos, end_pos, clade, n_snps) in enumerate(merged_tracts):
745770
cls = "ia" if clade == "ia" else "ib"
746771
lbl = ref1_label if clade == "ia" else ref2_label
772+
left_pct = start_pos / genome_length * 100 if genome_length else 0
773+
width_pct = max(0.3, (end_pos - start_pos + 1) / genome_length * 100) if genome_length else 2
747774
strip_segments += (
748-
f'<span class="strip-segment region-segment {cls}" title="{start_pos}{end_pos} {html_escape(lbl)} ({n_snps} SNPs)" style="flex: {n_snps} 1 0;"></span>'
775+
f'<span class="strip-segment region-segment {cls}" title="{start_pos}{end_pos} {html_escape(lbl)} ({n_snps} SNPs)" style="left:{left_pct:.3f}%; width:{width_pct:.3f}%;"></span>'
749776
)
750-
if j < len(merged_tracts) - 1:
751-
start_next = merged_tracts[j + 1][0]
752-
ca, cb = clade, merged_tracts[j + 1][2]
753-
lbl_a = ref1_label if ca == "ia" else ref2_label
754-
lbl_b = ref1_label if cb == "ia" else ref2_label
755-
bp_title = f"Breakpoint: {end_pos}{start_next} ({lbl_a}{lbl_b})"
756-
strip_segments += f'<span class="strip-segment breakpoint-marker" title="{html_escape(bp_title)}"></span>'
757777
n_tracts = len(merged_tracts)
758778
n_breakpoints = max(0, n_tracts - 1)
759779
section_cls = "rec-sites-section" + (" recombinant" if rec_call == "potential recombinant" else "")
@@ -767,8 +787,11 @@ def _ref_box(label: str, spec: Optional[str]) -> str:
767787
summary_text = "No recombination tracts (genome entirely one clade)"
768788
details_content = '<p class="threshold-note">No recombination detected; genome is entirely one clade.</p>'
769789
else:
770-
strip_min_w = max(400, (len(merged_tracts) + max(0, len(merged_tracts) - 1)) * 24)
771-
strip_display = f'<div class="strip-genome breakpoints-strip" style="min-width:{strip_min_w}px" role="img" aria-label="Predicted regions and breakpoints">{strip_segments}</div>'
790+
bp_ruler_html = _genome_ruler_html(genome_length, bp_strip_min_w)
791+
strip_display = (
792+
f'<div class="strip-genome breakpoints-strip" style="min-width:{bp_strip_min_w}px" role="img" aria-label="Predicted regions and breakpoints">{strip_segments}</div>'
793+
+ bp_ruler_html
794+
)
772795
summary_text = f"Show recombination tracts (Number of tracts: {n_tracts}, breakpoints: {n_breakpoints})"
773796
details_content = (
774797
f'<table class="rec-sites-table"><thead><tr><th>Tract #</th><th>Beginning of tract (bp)</th><th>End of tract (bp)</th><th>Clade</th></tr></thead><tbody>'
@@ -791,11 +814,11 @@ def _ref_box(label: str, spec: Optional[str]) -> str:
791814
strip_ref2_name = ref2_label if ref2_label not in ("ref1", "ref2") else (ref2_spec or ref2_label)
792815
breakpoints_section_html = (
793816
'<details class="collapsible-section diagnostic-strips-chart" open id="breakpointsStripsSection">'
794-
'<summary><h2>Recombination breakpoints per sample</h2></summary>'
817+
'<summary><h2>Recombination tracts and predicted breakpoints per sample</h2></summary>'
795818
'<div class="section-inner chart-section">'
796-
'<p class="threshold-note">Predicted recombination breakpoints within each genome. We show the beginning and end of each detected tract (first and last diagnostic SNP of that clade). The <strong>breakpoint lies in the region between</strong> the end of one tract and the start of the next; we cannot pinpoint its exact position because those regions have no diagnostic SNPs (genetically identical). Minimum consecutive diagnostic SNPs per tract: <strong>{min_consecutive}</strong>. <span id="breakpointsFilterCount" aria-live="polite"></span></p>'
797-
'<p class="threshold-note">{ref1} (blue), {ref2} (orange), breakpoint (red bar).</p>'
798-
'<div class="strip-legend"><span class="strip-legend-ia"></span> {ref1} &nbsp; <span class="strip-legend-ib"></span> {ref2} &nbsp; <span class="strip-legend-breakpoint"></span> breakpoint</div>'
819+
'<p class="threshold-note">Each coloured tract spans from the <strong>first to the last diagnostic SNP</strong> unambiguously derived from that clade. The predicted breakpoint lies somewhere in the <strong>uncoloured gap</strong> between adjacent tracts — its exact position cannot be determined because those intervening regions lack clade-informative diagnostic SNPs. Minimum consecutive diagnostic SNPs per tract: <strong>{min_consecutive}</strong>. <span id="breakpointsFilterCount" aria-live="polite"></span></p>'
820+
'<p class="threshold-note">{ref1} (blue), {ref2} (orange). Grey gaps = predicted breakpoint region (may be widened by ambiguous bases or poorly sequenced areas).</p>'
821+
'<div class="strip-legend"><span class="strip-legend-ia"></span> {ref1} &nbsp; <span class="strip-legend-ib"></span> {ref2} &nbsp; <span class="strip-legend-gap"></span> predicted breakpoint region (affected by ambiguous bases / poor coverage)</div>'
799822
'<div class="strip-strips-container" id="breakpointsStripScrollWrapper">'
800823
'<div id="breakpointsStripsContainer">'
801824
).format(ref1=html_escape(strip_ref1_name), ref2=html_escape(strip_ref2_name), min_consecutive=min_consecutive)
@@ -901,19 +924,27 @@ def _ref_box(label: str, spec: Optional[str]) -> str:
901924
.rec-sites-row {{ display: flex; align-items: center; gap: 12px; flex-wrap: nowrap; min-width: 0; width: 100%; max-width: 100%; }}
902925
.rec-sites-sample-id {{ font-size: 0.95em; font-weight: 600; color: #333; width: 300px; min-width: 300px; max-width: 300px; overflow: visible; white-space: normal; word-break: break-word; flex-shrink: 0; }}
903926
.strip-cell {{ flex: 1 0 0; min-width: 0; overflow-x: auto; overflow-y: hidden; border-radius: 4px; border: 1px solid #e9ecef; -webkit-overflow-scrolling: touch; }}
904-
.strip-genome {{ display: flex; flex-wrap: nowrap; height: 24px; min-width: 200px; border-radius: 4px; overflow: hidden; }}
905-
.strip-segment {{ flex: 1; min-width: 2px; transition: opacity 0.15s; }}
906-
.strip-segment:hover {{ opacity: 0.85; }}
927+
.strip-genome {{ position: relative; width: 100%; height: 24px; min-width: 200px; border-radius: 4px; overflow: hidden; background: #e9ecef; }}
928+
.strip-genome.breakpoints-strip {{ background: #ced4da; height: 32px; border-radius: 6px; box-shadow: inset 0 2px 6px rgba(0,0,0,0.13); }}
929+
.strip-genome.breakpoints-strip::after {{ content: ''; position: absolute; inset: 0; background: linear-gradient(to bottom, rgba(255,255,255,0.18) 0%, transparent 55%); pointer-events: none; z-index: 5; border-radius: inherit; }}
930+
.strip-segment {{ position: absolute; top: 0; height: 100%; width: 2px; transition: opacity 0.15s; }}
931+
.strip-segment:hover {{ opacity: 0.75; }}
907932
.strip-segment.ia {{ background: #4A90D9; }}
908933
.strip-segment.ib {{ background: #E89B3C; }}
909934
.strip-segment.other {{ background: #95a5a6; }}
935+
.strip-segment.region-segment {{ min-width: 4px; border-radius: 3px; }}
936+
.strip-segment.breakpoint-marker {{ width: 4px; background: #c0392b; transform: translateX(-50%); }}
937+
#breakpointsStripsContainer .rec-sites-section {{ background: linear-gradient(135deg, #fafbfc 0%, #f4f6f9 100%); border-left: 4px solid #dee2e6; box-shadow: 0 1px 4px rgba(0,0,0,0.05); transition: box-shadow 0.15s; }}
938+
#breakpointsStripsContainer .rec-sites-section.recombinant {{ border-left-color: #E89B3C; }}
939+
#breakpointsStripsContainer .rec-sites-section:hover {{ box-shadow: 0 3px 10px rgba(0,0,0,0.10); }}
910940
.strip-legend {{ display: flex; align-items: center; gap: 4px; flex-wrap: wrap; margin-bottom: 12px; font-size: 0.9em; color: #495057; }}
911941
.strip-legend-ia {{ display: inline-block; width: 14px; height: 14px; background: #4A90D9; border-radius: 2px; }}
912942
.strip-legend-ib {{ display: inline-block; width: 14px; height: 14px; background: #E89B3C; border-radius: 2px; }}
913943
.strip-legend-other {{ display: inline-block; width: 14px; height: 14px; background: #95a5a6; border-radius: 2px; }}
914-
.strip-legend-breakpoint {{ display: inline-block; width: 4px; height: 14px; background: #c0392b; border-radius: 1px; }}
915-
.strip-segment.breakpoint-marker {{ flex: none; width: 4px; min-width: 4px; background: #c0392b; }}
916-
.strip-genome.breakpoints-strip .strip-segment.region-segment {{ min-width: 8px; }}
944+
.strip-legend-gap {{ display: inline-block; width: 14px; height: 14px; background: #ced4da; border: 1px solid #adb5bd; border-radius: 2px; }}
945+
.strip-ruler {{ position: relative; width: 100%; height: 22px; min-width: 200px; margin-top: 3px; }}
946+
.ruler-tick {{ position: absolute; transform: translateX(-50%); font-size: 0.67em; font-weight: 500; color: #6c757d; white-space: nowrap; line-height: 1; padding-top: 6px; letter-spacing: 0.01em; }}
947+
.ruler-tick::before {{ content: ''; display: block; position: absolute; top: 0; left: 50%; transform: translateX(-50%); width: 1px; height: 5px; background: #adb5bd; }}
917948
.rec-sites-details {{ margin-top: 10px; font-size: 0.9em; }}
918949
.rec-sites-details summary {{ cursor: pointer; color: #667eea; font-weight: 500; }}
919950
.rec-sites-table {{ margin-top: 8px; border-collapse: collapse; font-size: 0.9em; max-height: 200px; overflow: auto; }}
@@ -1120,12 +1151,13 @@ def _ref_box(label: str, spec: Optional[str]) -> str:
11201151
}});
11211152
}});
11221153
1154+
function truncId(id) {{ return id && id.length > 20 ? id.slice(0, 20) + "\u2026" : (id || ""); }}
11231155
var chartBar = null;
11241156
function updateChart() {{
11251157
var visibleRows = rows.filter(function(r) {{ return !r.classList.contains("hidden"); }});
11261158
var labels = visibleRows.map(function(r) {{
11271159
var ri = parseInt(r.getAttribute("data-row"), 10);
1128-
return (data[ri] && data[ri].id) ? data[ri].id : "";
1160+
return truncId((data[ri] && data[ri].id) ? data[ri].id : "");
11291161
}});
11301162
var pct1 = visibleRows.map(function(r) {{
11311163
var ri = parseInt(r.getAttribute("data-row"), 10);
@@ -1171,7 +1203,7 @@ def _ref_box(label: str, spec: Optional[str]) -> str:
11711203
}}
11721204
}},
11731205
scales: {{
1174-
x: {{ title: {{ display: true, text: "Accession" }}, ticks: {{ maxRotation: 45, minRotation: 45, autoSkip: false, font: {{ size: 11 }} }} }},
1206+
x: {{ title: {{ display: true, text: "Sample ID" }}, ticks: {{ maxRotation: 45, minRotation: 45, autoSkip: false, font: {{ size: 11 }} }} }},
11751207
y: {{ title: {{ display: true, text: "Percentage (%)" }}, min: 0, max: 100, ticks: {{ stepSize: 20 }} }}
11761208
}},
11771209
plugins: {{ legend: {{ display: true, position: "top" }} }}

0 commit comments

Comments
 (0)