"""Off-target analysis for siRNA design.
This module provides comprehensive off-target analysis functionality for siRNA design,
including both miRNA seed match analysis and transcriptome off-target detection.
Uses BWA-MEM2 for transcriptome alignments and supports in-process miRNA seed scanning.
Optimized for both standalone use and parallelized Nextflow workflows.
"""
import importlib
import json
import os
import re
import shutil
import statistics
import subprocess # nosec B404
import tempfile
from collections import defaultdict
from dataclasses import dataclass
from enum import Enum
from pathlib import Path
from typing import Any, Protocol
import pandas as pd
from sirnaforge.data.base import FastaUtils
from sirnaforge.data.mirna_manager import MiRNADatabaseManager
from sirnaforge.models.off_target import (
AggregatedMiRNASummary,
AggregatedOffTargetSummary,
AlignmentStrand,
AnalysisMode,
AnalysisSummary,
MiRNAHit,
MiRNASummary,
OffTargetHit,
)
from sirnaforge.models.schemas import GenomeAlignmentSchema, MiRNAAlignmentSchema
from sirnaforge.models.sirna import SiRNACandidate
from sirnaforge.utils.logging_utils import get_logger
from sirnaforge.utils.species import human_vs_other_totals
logger = get_logger(__name__)
_NUCLEOTIDE_ALPHABET = ("A", "C", "G", "T")
_CIGAR_OP_PATTERN = re.compile(r"(\d+)([MIDNSHP=X])")
def _mirna_max_hits() -> int | None:
"""Resolve the miRNA seed-hit cap (per species, batched across candidates).
Exhaustive by default (``None`` = report every seed match) — a truncated cap
silently loses hits and biases downstream counts. Set
``SIRNAFORGE_MIRNA_MAX_HITS`` to a positive int to impose a cap if desired.
"""
raw = os.getenv("SIRNAFORGE_MIRNA_MAX_HITS")
if raw:
try:
value = int(raw)
if value > 0:
return value
logger.warning("SIRNAFORGE_MIRNA_MAX_HITS must be positive; ignoring (exhaustive).")
except ValueError:
logger.warning("Invalid SIRNAFORGE_MIRNA_MAX_HITS=%r; ignoring (exhaustive).", raw)
return None
[docs]
class MiRNASeedBackend(str, Enum):
"""Internal backend options for miRNA seed scanning."""
BWA = "bwa"
EXHAUSTIVE_PYTHON = "exhaustive_python"
PYAHOCORASICK = "pyahocorasick"
@dataclass(frozen=True, slots=True)
class _PreparedSeedQuery:
"""Normalized siRNA query plus extracted seed sequence."""
qname: str
qseq: str
seed_qseq: str
@dataclass(frozen=True, slots=True)
class _MiRNASeedMatch:
"""Semantic seed-match hit used by in-process backends."""
qname: str
qseq: str
seed_qseq: str
mirna_id: str
coord: int
mismatch_positions: tuple[int, ...]
@dataclass(frozen=True, slots=True)
class _NormalizedMiRNASeedHit:
"""Backend-independent miRNA seed hit contract used by adapters and parity checks."""
qname: str
qseq: str
mirna_id: str
coord: int
strand: str
cigar: str
mapq: int
as_score: int | None
nm: int
mismatch_positions: tuple[int, ...]
seed_mismatches: int
offtarget_score: float
def to_alignment_row(self, *, species: str, database: str) -> dict[str, Any]:
"""Project the normalized hit onto the stable MiRNAAlignmentSchema row contract."""
return {
"qname": self.qname,
"qseq": self.qseq,
"species": species,
"database": database,
"mirna_id": self.mirna_id,
"coord": self.coord,
"strand": self.strand,
"cigar": self.cigar,
"mapq": self.mapq,
"as_score": self.as_score,
"nm": self.nm,
"seed_mismatches": self.seed_mismatches,
"offtarget_score": self.offtarget_score,
}
def semantic_identity(self) -> tuple[str, str, int, int, int, float]:
"""Return the backend-independent semantic identity used for parity checks."""
return (
self.qname,
self.mirna_id,
self.coord,
self.nm,
self.seed_mismatches,
self.offtarget_score,
)
class _MiRNASeedScanner(Protocol):
"""Internal contract for miRNA seed scanners."""
backend: MiRNASeedBackend
def scan(
self,
queries: list[_PreparedSeedQuery],
mirna_sequences: dict[str, str],
*,
max_mismatches: int,
) -> list[_MiRNASeedMatch]:
"""Return semantic seed hits across one miRNA FASTA dictionary."""
...
class _MiRNASeedBackendUnavailableError(RuntimeError):
"""Raised when an in-process backend dependency is unavailable."""
def _compute_species_counts(df: pd.DataFrame) -> dict[str, int]:
"""Build a frequency map for the 'species' column in a DataFrame."""
if df.empty or "species" not in df.columns:
return {}
counts: dict[str, int] = {}
for value in df["species"].tolist():
label = "unknown"
if value is not None and not pd.isna(value):
label = str(value)
counts[label] = counts.get(label, 0) + 1
return counts
def _normalize_nucleotide_sequence(sequence: str) -> str:
"""Normalize RNA/DNA sequences into the comparison alphabet used by seed scanners."""
return sequence.upper().replace("U", "T")
def _optional_int(value: Any) -> int | None:
"""Convert loose backend metadata to int while preserving explicit missing sentinels."""
if value in (None, "", "NA"):
return None
return int(value)
def _prepare_seed_queries(
sequences: dict[str, str],
*,
seed_start: int,
seed_end: int,
) -> list[_PreparedSeedQuery]:
"""Extract normalized seed queries from full siRNA guide sequences."""
prepared: list[_PreparedSeedQuery] = []
for name, sequence in sequences.items():
normalized = _normalize_nucleotide_sequence(sequence)
if len(normalized) >= seed_end:
seed_qseq = normalized[seed_start - 1 : seed_end]
else:
logger.warning(f"Sequence {name} too short for seed extraction: {sequence}")
seed_qseq = normalized
prepared.append(_PreparedSeedQuery(qname=name, qseq=normalized, seed_qseq=seed_qseq))
return prepared
def _mismatch_positions(query: str, window: str) -> tuple[int, ...]:
"""Return 1-based mismatch positions for two equal-length sequences."""
return tuple(index + 1 for index, (left, right) in enumerate(zip(query, window, strict=True)) if left != right)
def _seed_window_positions_to_guide(
window_positions: tuple[int, ...] | list[int],
*,
seed_start: int,
) -> tuple[int, ...]:
"""Map 1-based positions within the extracted seed window onto 1-based guide positions.
The in-process seed scanners align only the guide seed window (guide positions
``seed_start``..``seed_end``) against the miRNA, so ``_mismatch_positions`` returns
positions relative to that window (1..len(seed)). Downstream seed/position logic
(``seed_mismatches`` counting and ``_calculate_seed_offtarget_score``) is expressed
in guide coordinates, so window position ``w`` must be shifted by ``seed_start - 1``.
"""
offset = seed_start - 1
return tuple(offset + position for position in window_positions)
@dataclass(frozen=True)
class _QueryLayout:
"""Where every base of one SAM record's read sits, in full-read coordinates.
Read coordinates include hard-clipped bases so offsets stay comparable to the original
query the aligner was handed. Three separate facts have to come out of the CIGAR before
an MD tag means anything in read coordinates:
* MD offsets count only reference-aligned read bases (``M``/``=``/``X``), so an insertion
or a leading clip desynchronises them from read positions -- hence ``md_to_read``.
* clipped and inserted read bases never appear in the MD tag at all, yet they failed to
pair with the target, so they carry off-target risk -- hence ``clipped``/``inserted``.
* ``D``/``N`` skip *reference* bases, so there is no read position to blame; the penalty
is anchored on the read base immediately 5' of the gap -- hence ``deletions``.
"""
read_length: int
md_to_read: tuple[int, ...]
"""1-based MD offset -> 1-based read position (index 0 is an unused placeholder)."""
clipped: tuple[int, ...]
inserted: tuple[int, ...]
deletions: tuple[tuple[int, int], ...]
"""``(anchor read position, skipped reference bases)`` for each ``D``/``N`` op."""
@property
def unpaired(self) -> tuple[int, ...]:
"""Read positions with no target base opposite them (clipped or inserted)."""
return self.clipped + self.inserted
def _parse_cigar_query_layout(cigar: str, *, read_length: int = 0) -> _QueryLayout:
"""Map one CIGAR onto read coordinates, keeping clips, insertions and gaps separate.
``read_length`` is the length of the query the aligner was handed; it is used as the
fallback layout when the record carries no CIGAR (``*``), in which case the whole read
is assumed to have aligned ungapped and MD offsets are already read positions.
"""
md_to_read: list[int] = [0] # 1-based indexing; slot 0 is never read
clipped: list[int] = []
inserted: list[int] = []
deletions: list[tuple[int, int]] = []
if not cigar or cigar == "*":
md_to_read.extend(range(1, read_length + 1))
return _QueryLayout(read_length, tuple(md_to_read), (), (), ())
position = 1
for length_str, op in _CIGAR_OP_PATTERN.findall(cigar):
length = int(length_str)
if op in "M=X":
md_to_read.extend(range(position, position + length))
position += length
elif op in "SH":
clipped.extend(range(position, position + length))
position += length
elif op == "I":
inserted.extend(range(position, position + length))
position += length
elif op in "DN":
deletions.append((max(1, position - 1), length))
# P (padding) consumes neither read nor reference bases.
return _QueryLayout(position - 1, tuple(md_to_read), tuple(clipped), tuple(inserted), tuple(deletions))
def _calculate_seed_offtarget_score(
mismatch_positions: list[int] | tuple[int, ...],
*,
seed_start: int,
seed_end: int,
) -> float:
"""Calculate the same mismatch-weighted penalty used by BWA-derived seed hits."""
num_mismatches = len(mismatch_positions)
if num_mismatches == 0:
return 0.0
base_score = num_mismatches * 10.0
position_penalty = 0.0
for pos in mismatch_positions:
if seed_start <= pos <= seed_end:
position_penalty += 5.0
elif pos <= 10:
position_penalty += 3.0
else:
position_penalty += 1.0
continuous_bonus = 0.0
if num_mismatches >= 2:
sorted_positions = sorted(mismatch_positions)
continuous_count = 0
for index in range(len(sorted_positions) - 1):
if sorted_positions[index + 1] - sorted_positions[index] == 1:
continuous_count += 1
continuous_bonus = continuous_count * 2.0
return base_score + position_penalty + continuous_bonus
def _alignment_score(query_length: int, mismatch_count: int) -> int:
"""Approximate ungapped seed alignment score used for audit metadata."""
return max(0, (query_length * 2) - (mismatch_count * 3))
[docs]
def mirna_seed_hit_identity(
hit: dict[str, Any],
*,
coord_is_one_based: bool = False,
) -> tuple[str, str, int, int, int, float]:
"""Return the semantic identity tuple used for backend comparison.
This normalizes the current in-process seed-scan rows and the existing
BWA-derived rows onto the same comparison contract.
"""
return normalize_mirna_seed_hit(hit, coord_is_one_based=coord_is_one_based).semantic_identity()
[docs]
def normalize_mirna_seed_hit(
hit: dict[str, Any],
*,
coord_is_one_based: bool = False,
) -> _NormalizedMiRNASeedHit:
"""Normalize backend-specific miRNA seed hit metadata onto one adapter/parity contract."""
mirna_id = str(hit.get("mirna_id") or hit.get("rname") or "")
coord_value = hit.get("coord", 0)
if isinstance(coord_value, str) and ":" in coord_value:
coord = int(str(coord_value).rsplit(":", maxsplit=1)[-1])
else:
coord = int(coord_value)
if coord_is_one_based:
coord -= 1
mismatch_positions_value = hit.get("mismatch_positions", ())
if mismatch_positions_value in (None, "", "NA"):
mismatch_positions: tuple[int, ...] = ()
elif isinstance(mismatch_positions_value, str):
stripped = mismatch_positions_value.strip().strip("[]()")
mismatch_positions = tuple(int(part.strip()) for part in stripped.split(",") if part.strip())
else:
mismatch_positions = tuple(int(position) for position in mismatch_positions_value)
as_score = _optional_int(hit.get("as_score"))
nm = _optional_int(hit.get("nm"))
if nm is None:
nm = len(mismatch_positions)
seed_mismatches = _optional_int(hit.get("seed_mismatches"))
if seed_mismatches is None:
seed_mismatches = nm
return _NormalizedMiRNASeedHit(
qname=str(hit["qname"]),
qseq=_normalize_nucleotide_sequence(str(hit.get("qseq", ""))),
mirna_id=mirna_id,
coord=coord,
strand=str(hit.get("strand", "+")),
cigar=str(hit.get("cigar", "")),
mapq=int(hit.get("mapq", 255)),
as_score=as_score,
nm=nm,
mismatch_positions=mismatch_positions,
seed_mismatches=seed_mismatches,
offtarget_score=float(hit.get("offtarget_score", 0.0)),
)
def _build_mirna_alignment_frame(
hits: list[dict[str, Any]],
*,
species: str,
database: str,
coord_is_one_based: bool = False,
) -> pd.DataFrame:
"""Adapt backend-specific hit rows into the stable miRNA alignment table contract."""
schema_columns = list(MiRNAAlignmentSchema.to_schema().columns.keys())
rows = [
normalize_mirna_seed_hit(hit, coord_is_one_based=coord_is_one_based).to_alignment_row(
species=species,
database=database,
)
for hit in hits
]
df: pd.DataFrame = pd.DataFrame(rows, columns=schema_columns)
return df
def _normalize_bwa_mirna_seed_hits(
hits: list[dict[str, Any]],
*,
sequences: dict[str, str],
mirna_sequences: dict[str, str],
seed_start: int,
seed_end: int,
max_mismatches: int,
) -> list[dict[str, Any]]:
"""Project BWA seed alignments onto the backend-independent semantic seed-hit contract."""
prepared_queries = {
query.qname: query for query in _prepare_seed_queries(sequences, seed_start=seed_start, seed_end=seed_end)
}
normalized_hits: dict[tuple[str, str, int, tuple[int, ...]], dict[str, Any]] = {}
for hit in hits:
qname = str(hit["qname"])
prepared_query = prepared_queries.get(qname)
mirna_id = str(hit.get("rname", ""))
mirna_sequence = mirna_sequences.get(mirna_id)
if prepared_query is None or mirna_sequence is None:
continue
coord_value = str(hit.get("coord", ""))
start = int(coord_value.rsplit(":", maxsplit=1)[-1]) - 1 if ":" in coord_value else int(coord_value) - 1
if start < 0:
continue
window_end = start + len(prepared_query.seed_qseq)
if window_end > len(mirna_sequence):
continue
target_window = _normalize_nucleotide_sequence(mirna_sequence[start:window_end])
window_positions = _mismatch_positions(prepared_query.seed_qseq, target_window)
if len(window_positions) > max_mismatches:
continue
# Positions are window-relative (1..len(seed)); shift onto guide coordinates so
# seed classification and scoring match the guide-relative seed_start/seed_end frame.
mismatch_positions = _seed_window_positions_to_guide(window_positions, seed_start=seed_start)
seed_mismatches = sum(1 for pos in mismatch_positions if seed_start <= pos <= seed_end)
normalized_hits[(qname, mirna_id, start, mismatch_positions)] = {
"qname": qname,
"qseq": prepared_query.qseq,
"rname": mirna_id,
"coord": start,
"strand": "+",
"cigar": f"{len(prepared_query.seed_qseq)}M",
"mapq": 255,
"as_score": _alignment_score(len(prepared_query.seed_qseq), len(mismatch_positions)),
"nm": len(mismatch_positions),
"mismatch_positions": list(mismatch_positions),
"seed_mismatches": seed_mismatches,
"offtarget_score": _calculate_seed_offtarget_score(
mismatch_positions,
seed_start=seed_start,
seed_end=seed_end,
),
}
results = list(normalized_hits.values())
results.sort(
key=lambda row: (
row["offtarget_score"],
-int(row.get("as_score") or 0),
row["qname"],
row["rname"],
row["coord"],
)
)
return results
def _candidate_patterns(query: str, max_mismatches: int) -> dict[str, tuple[int, ...]]:
"""Generate exact candidate patterns up to a maximum mismatch count."""
patterns: dict[str, tuple[int, ...]] = {query: ()}
query_chars = list(query)
def walk(index: int, mismatch_positions: tuple[int, ...], current: list[str]) -> None:
if len(mismatch_positions) > max_mismatches:
return
if index == len(query_chars):
candidate = "".join(current)
existing = patterns.get(candidate)
if existing is None or len(mismatch_positions) < len(existing):
patterns[candidate] = mismatch_positions
return
original = query_chars[index]
current.append(original)
walk(index + 1, mismatch_positions, current)
current.pop()
if len(mismatch_positions) == max_mismatches:
return
for base in _NUCLEOTIDE_ALPHABET:
if base == original:
continue
current.append(base)
walk(index + 1, mismatch_positions + (index + 1,), current)
current.pop()
walk(0, (), [])
return patterns
class _ExhaustivePythonMiRNASeedScanner:
"""Baseline seed scanner used as the correctness oracle in tests."""
backend = MiRNASeedBackend.EXHAUSTIVE_PYTHON
def scan(
self,
queries: list[_PreparedSeedQuery],
mirna_sequences: dict[str, str],
*,
max_mismatches: int,
) -> list[_MiRNASeedMatch]:
matches: list[_MiRNASeedMatch] = []
normalized_mirnas = {mirna_id: _normalize_nucleotide_sequence(seq) for mirna_id, seq in mirna_sequences.items()}
for query in queries:
query_length = len(query.seed_qseq)
if query_length == 0:
continue
for mirna_id, mirna_seq in normalized_mirnas.items():
if len(mirna_seq) < query_length:
continue
for start in range(len(mirna_seq) - query_length + 1):
window = mirna_seq[start : start + query_length]
mismatch_positions = _mismatch_positions(query.seed_qseq, window)
if len(mismatch_positions) > max_mismatches:
continue
matches.append(
_MiRNASeedMatch(
qname=query.qname,
qseq=query.qseq,
seed_qseq=query.seed_qseq,
mirna_id=mirna_id,
coord=start,
mismatch_positions=mismatch_positions,
)
)
return matches
class _PyAhoCorasickMiRNASeedScanner:
"""pyahocorasick-based scanner using pre-expanded mismatch patterns."""
backend = MiRNASeedBackend.PYAHOCORASICK
@staticmethod
def _load_module() -> Any:
try:
return importlib.import_module("ahocorasick")
except ImportError as exc:
raise _MiRNASeedBackendUnavailableError(
"miRNA seed backend 'pyahocorasick' is unavailable because the installed environment is missing "
"the required 'pyahocorasick' package"
) from exc
def scan(
self,
queries: list[_PreparedSeedQuery],
mirna_sequences: dict[str, str],
*,
max_mismatches: int,
) -> list[_MiRNASeedMatch]:
ahocorasick = self._load_module()
automaton = ahocorasick.Automaton()
payloads_by_pattern: dict[str, list[_PreparedSeedQuery]] = defaultdict(list)
for query in queries:
if not query.seed_qseq:
continue
for pattern in _candidate_patterns(query.seed_qseq, max_mismatches):
payloads_by_pattern[pattern].append(query)
for pattern, pattern_queries in payloads_by_pattern.items():
automaton.add_word(pattern, (pattern, tuple(pattern_queries)))
automaton.make_automaton()
matches: list[_MiRNASeedMatch] = []
normalized_mirnas = {mirna_id: _normalize_nucleotide_sequence(seq) for mirna_id, seq in mirna_sequences.items()}
for mirna_id, mirna_seq in normalized_mirnas.items():
for end_index, (pattern, pattern_queries) in automaton.iter(mirna_seq):
start = end_index - len(pattern) + 1
window = mirna_seq[start : end_index + 1]
for query in pattern_queries:
mismatch_positions = _mismatch_positions(query.seed_qseq, window)
if len(mismatch_positions) > max_mismatches:
continue
matches.append(
_MiRNASeedMatch(
qname=query.qname,
qseq=query.qseq,
seed_qseq=query.seed_qseq,
mirna_id=mirna_id,
coord=start,
mismatch_positions=mismatch_positions,
)
)
return matches
def _get_mirna_seed_scanner(backend: MiRNASeedBackend) -> _MiRNASeedScanner:
"""Resolve one in-process miRNA seed backend implementation."""
if backend == MiRNASeedBackend.EXHAUSTIVE_PYTHON:
return _ExhaustivePythonMiRNASeedScanner()
if backend == MiRNASeedBackend.PYAHOCORASICK:
return _PyAhoCorasickMiRNASeedScanner()
raise ValueError(f"Unsupported in-process miRNA seed backend: {backend.value}")
[docs]
def scan_mirna_seed_matches(
sequences: dict[str, str],
mirna_sequences: dict[str, str],
*,
backend: MiRNASeedBackend | str = MiRNASeedBackend.PYAHOCORASICK,
seed_start: int = 2,
seed_end: int = 8,
max_mismatches: int = 2,
max_hits: int | None = None,
) -> list[dict[str, Any]]:
"""Scan miRNA FASTA records using an in-process seed-scanning backend.
The returned rows preserve the existing internal hit shape used by miRNA analysis,
so callers can continue adapting them to DataFrame-based outputs.
"""
resolved_backend = MiRNASeedBackend(backend)
if resolved_backend == MiRNASeedBackend.BWA:
raise ValueError("scan_mirna_seed_matches only supports in-process backends")
queries = _prepare_seed_queries(sequences, seed_start=seed_start, seed_end=seed_end)
scanner = _get_mirna_seed_scanner(resolved_backend)
raw_matches = scanner.scan(queries, mirna_sequences, max_mismatches=max_mismatches)
deduped: dict[tuple[str, str, int, tuple[int, ...]], dict[str, Any]] = {}
for match in raw_matches:
# ``_mismatch_positions`` reports positions within the extracted seed window
# (1..len(seed)); shift them onto guide coordinates so seed classification and
# position-weighted scoring are computed in the same frame as the seed_start/
# seed_end thresholds.
mismatch_positions = _seed_window_positions_to_guide(match.mismatch_positions, seed_start=seed_start)
nm = len(mismatch_positions)
seed_mismatches = sum(1 for pos in mismatch_positions if seed_start <= pos <= seed_end)
key = (match.qname, match.mirna_id, match.coord, mismatch_positions)
deduped[key] = {
"qname": match.qname,
"qseq": match.qseq,
"rname": match.mirna_id,
"coord": match.coord,
"strand": "+",
"cigar": f"{len(match.seed_qseq)}M",
"mapq": 255,
"as_score": _alignment_score(len(match.seed_qseq), nm),
"nm": nm,
"mismatch_positions": list(mismatch_positions),
"seed_mismatches": seed_mismatches,
"offtarget_score": _calculate_seed_offtarget_score(
mismatch_positions,
seed_start=seed_start,
seed_end=seed_end,
),
}
results = list(deduped.values())
results.sort(
key=lambda row: (
row["offtarget_score"],
-int(row.get("as_score") or 0),
row["qname"],
row["rname"],
row["coord"],
)
)
return results if max_hits is None else results[:max_hits]
# =============================================================================
# Core Analyzer Classes
# =============================================================================
def _get_executable_path(tool_name: str) -> str | None:
"""Get the full path to an executable, ensuring it exists."""
path = shutil.which(tool_name)
if path is None:
logger.warning(f"Tool '{tool_name}' not found in PATH")
return path
def _validate_command_args(cmd: list[str]) -> None:
"""Validate command arguments for subprocess execution."""
if not cmd:
raise ValueError("Command list cannot be empty")
executable = cmd[0]
if not executable:
raise ValueError("Executable path cannot be empty")
# Ensure we have an absolute path to the executable
if not Path(executable).is_absolute():
raise ValueError(f"Executable must be an absolute path: {executable}")
# =============================================================================
# Core Analyzer Classes
# =============================================================================
[docs]
class BwaAnalyzer:
"""BWA-MEM2 based analyzer for both transcriptome and miRNA seed off-target search."""
[docs]
def __init__(
self,
index_prefix: str | Path,
mode: str = "transcriptome", # "transcriptome" or "mirna_seed"
seed_length: int = 12,
min_score: int = 15,
max_hits: int | None = None,
seed_start: int = 2,
seed_end: int = 8,
):
"""Initialize BWA-MEM2 analyzer.
Args:
index_prefix: Path to BWA index
mode: Analysis mode - "transcriptome" for long targets, "mirna_seed" for short targets
seed_length: BWA seed length parameter
min_score: Minimum alignment score
max_hits: Maximum hits to return (``None`` = no limit / exhaustive)
seed_start: Seed region start (1-based)
seed_end: Seed region end (1-based)
"""
self.index_prefix = str(index_prefix)
self.mode = mode
self.seed_length = seed_length
self.min_score = min_score
self.max_hits = max_hits
self.seed_start = seed_start
self.seed_end = seed_end
# Configure parameters based on mode
if mode == "mirna_seed":
# For miRNA seed analysis: short query (6-8bp) vs short target (~22bp)
# Need very permissive parameters for ultra-short sequences
self.seed_length = min(seed_length, 6) # Max 6bp seed for 6-8bp queries
self.min_score = 6 # Very low threshold - allow imperfect matches
elif mode == "transcriptome":
# For transcriptome analysis: short query vs long target
self.seed_length = seed_length # Use provided seed length
self.min_score = min_score # Use provided min score
else:
raise ValueError(f"Unknown mode: {mode}. Use 'transcriptome' or 'mirna_seed'")
[docs]
def analyze_sequences(self, sequences: dict[str, str]) -> list[dict[str, Any]]:
"""Run BWA-MEM2 analysis on sequences.
Args:
sequences: Dictionary of sequence name -> sequence
Returns:
List of alignment dictionaries
"""
# Prepare sequences based on mode
analysis_sequences = self._prepare_sequences_for_analysis(sequences)
results = []
temp_fasta_path = create_temp_fasta(analysis_sequences)
try:
# Get absolute path to bwa-mem2 executable
bwa_path = _get_executable_path("bwa-mem2")
if not bwa_path:
raise FileNotFoundError("BWA-MEM2 executable not found in PATH")
# Configure BWA parameters based on mode
cmd = self._build_bwa_command(bwa_path, temp_fasta_path)
_validate_command_args(cmd)
logger.info(f"Running BWA-MEM2 ({self.mode} mode): {' '.join(cmd)}")
result = subprocess.run(cmd, capture_output=True, text=True, timeout=None, check=True) # nosec B603
results = self._parse_sam_output(result.stdout, sequences)
results = self._filter_and_rank(results)
logger.info(f"BWA-MEM2 analysis completed: {len(results)} hits found")
except subprocess.CalledProcessError as e:
logger.error(f"BWA-MEM2 failed: {e.stderr}")
except subprocess.TimeoutExpired:
logger.error("BWA-MEM2 timed out")
finally:
Path(temp_fasta_path).unlink(missing_ok=True)
return results if self.max_hits is None else results[: self.max_hits]
def _prepare_sequences_for_analysis(self, sequences: dict[str, str]) -> dict[str, str]:
"""Prepare sequences for analysis based on mode."""
if self.mode == "mirna_seed":
prepared = _prepare_seed_queries(sequences, seed_start=self.seed_start, seed_end=self.seed_end)
return {query.qname: query.seed_qseq for query in prepared}
# For transcriptome mode, use full sequences
return {name: _normalize_nucleotide_sequence(seq) for name, seq in sequences.items()}
def _build_bwa_command(self, bwa_path: str, temp_fasta_path: str) -> list[str]:
"""Build BWA command based on analysis mode."""
base_cmd = [
bwa_path,
"mem",
"-a", # Output all alignments
"-v",
"1", # Verbosity level
]
if self.mode == "mirna_seed":
# For miRNA seed analysis: ultra-permissive parameters for 6-8bp vs ~22bp
cmd = base_cmd + [
"-k",
str(self.seed_length), # Seed length (max 6bp)
"-T",
str(self.min_score), # Minimum score (6)
"-w",
"2", # Narrow band width for short sequences
"-A",
"2", # Higher matching score to reward matches
"-B",
"1", # Low mismatch penalty
"-O",
"1,1", # Low gap open penalties
"-E",
"1,1", # Low gap extension penalties
"-L",
"8,8", # Clipping penalty for ultra-short reads
self.index_prefix,
temp_fasta_path,
]
elif self.mode == "transcriptome":
# For transcriptome analysis: standard parameters
cmd = base_cmd + [
"-k",
str(self.seed_length), # Seed length
"-T",
str(self.min_score), # Minimum score
"-w",
"100", # Band width (larger for long targets)
self.index_prefix,
temp_fasta_path,
]
else:
raise ValueError(f"Unknown mode: {self.mode}")
return cmd
def _query_frames(self, original_sequences: dict[str, str]) -> dict[str, tuple[int, int]]:
"""Map each query name onto ``(query length, guide-coordinate offset)``.
The aligner is handed a different query per mode (see
``_prepare_sequences_for_analysis``): the whole normalized guide in ``transcriptome``
mode, but only the extracted seed window in ``mirna_seed`` mode. Read positions in the
latter are therefore *window* positions and must be shifted by ``seed_start - 1`` before
they can be compared against the guide-relative ``seed_start``/``seed_end`` thresholds --
the same frame correction ``_normalize_bwa_mirna_seed_hits`` and
``scan_mirna_seed_matches`` already apply to their own rows. The lengths also give the
true read length for records whose CIGAR is absent (``*``).
"""
if self.mode == "mirna_seed":
return {
query.qname: (
len(query.seed_qseq),
# _prepare_seed_queries falls back to the whole guide when the guide is
# shorter than seed_end; read positions are then already guide positions.
self.seed_start - 1 if len(query.qseq) >= self.seed_end else 0,
)
for query in _prepare_seed_queries(
original_sequences, seed_start=self.seed_start, seed_end=self.seed_end
)
}
return {name: (len(_normalize_nucleotide_sequence(seq)), 0) for name, seq in original_sequences.items()}
def _guide_frame_mismatches(
self,
*,
flag: int,
cigar: str,
md_tag: str,
edit_distance: int,
query_len: int,
guide_offset: int,
) -> tuple[list[int], int, float]:
"""Return ``(mismatch_positions, nm, offtarget_score)`` in guide coordinates.
``mismatch_positions`` are the 1-based *guide* positions that failed to pair with the
target: MD substitutions plus every clipped or inserted guide base. Reference bases
skipped by a ``D``/``N`` op have no guide position of their own, so they only enter the
score, anchored on the guide base flanking the bulge. ``nm`` counts all of those
mismatch-equivalents (see ``OffTargetHit.nm``), which is why it can exceed the
aligner's ``NM`` tag.
"""
layout = _parse_cigar_query_layout(cigar, read_length=query_len)
read_length = layout.read_length or query_len
if query_len and layout.read_length and layout.read_length != query_len:
logger.warning(
f"CIGAR '{cigar}' covers {layout.read_length} read bases but the aligner was handed "
f"{query_len}; using the CIGAR length for guide-frame mapping"
)
md_offsets = self._parse_md_tag(md_tag)
substitutions = [layout.md_to_read[offset] for offset in md_offsets if 0 < offset < len(layout.md_to_read)]
if len(substitutions) != len(md_offsets):
logger.warning(
f"MD tag '{md_tag}' does not fit CIGAR '{cigar}': "
f"{len(md_offsets) - len(substitutions)} mismatch position(s) dropped"
)
# BWA stores flag&16 records reverse-complemented, so SEQ/MD/CIGAR are in the revcomp
# frame while the seed window is defined in guide coordinates: read position r is guide
# position read_length + 1 - r. design.py builds every guide as
# reverse_complement(target_seq), which makes minus-strand the common case, not an edge.
if flag & 16 and not read_length:
logger.warning(
f"Minus-strand record with no CIGAR and unknown query length (MD '{md_tag}'): mismatch "
"positions stay in the read frame, so seed classification for this hit is unreliable"
)
def to_guide(read_positions: list[int]) -> list[int]:
mirrored = (
[read_length + 1 - pos for pos in read_positions] if (flag & 16 and read_length) else read_positions
)
return sorted(guide_offset + pos for pos in mirrored)
unpaired = list(layout.unpaired)
scoring_reads = substitutions + unpaired
for anchor, gap_length in layout.deletions:
scoring_reads.extend([anchor] * gap_length)
mismatch_positions = to_guide(substitutions + unpaired)
# Never fall below what the aligner itself reported, so a missing or garbled MD tag can
# only ever over-penalise a hit.
nm = max(len(scoring_reads), edit_distance + len(layout.clipped))
offtarget_score = self._calculate_offtarget_score(to_guide(scoring_reads))
if nm and not offtarget_score:
# Invariant: only a full-length exact match may score 0.0 (= highest risk), because
# _filter_and_rank sorts ASCENDING and max_hits truncation keeps the head of that
# list. When differences are known to exist but could not be placed on the guide,
# charge the position-agnostic base term instead of letting the hit sort ahead of
# genuine perfect matches.
offtarget_score = nm * 10.0
return mismatch_positions, nm, offtarget_score
def _parse_sam_output(self, sam_output: str, original_sequences: dict[str, str]) -> list[dict[str, Any]]:
"""Parse SAM output from BWA-MEM2."""
results = []
query_frames = self._query_frames(original_sequences)
for line in sam_output.splitlines():
if line.startswith("@") or not line.strip():
continue
parts = line.split("\t")
if len(parts) < 11:
continue
qname = parts[0]
flag = int(parts[1])
rname = parts[2]
pos = int(parts[3])
mapq = int(parts[4]) if parts[4] != "*" else 0
cigar = parts[5]
if flag & 4: # Skip unmapped
continue
strand = "-" if (flag & 16) else "+"
coord = f"{rname}:{pos}"
# Parse optional tags
tags = {}
for tag in parts[11:]:
if ":" in tag:
tag_parts = tag.split(":", 2)
if len(tag_parts) == 3:
tags[tag_parts[0]] = tag_parts[2]
as_score = int(tags.get("AS", 0)) if "AS" in tags else None
query_len, guide_offset = query_frames.get(qname, (0, 0))
if not query_len and parts[9] != "*":
# SEQ is the read as the aligner stored it, so its length is the read length for
# everything but hard-clipped records -- a last resort for an unrecognised qname.
query_len = len(parts[9])
mismatch_positions, nm, offtarget_score = self._guide_frame_mismatches(
flag=flag,
cigar=cigar,
md_tag=tags.get("MD", ""),
edit_distance=int(tags.get("NM", 0)),
query_len=query_len,
guide_offset=guide_offset,
)
seed_mismatches = sum(1 for pos in mismatch_positions if self.seed_start <= pos <= self.seed_end)
result = {
"qname": qname,
"qseq": original_sequences.get(qname, ""), # Use original sequence
"rname": rname,
"coord": coord,
"strand": strand,
"cigar": cigar,
"mapq": mapq,
"as_score": as_score,
"nm": nm,
"mismatch_positions": mismatch_positions,
"seed_mismatches": seed_mismatches,
"offtarget_score": offtarget_score,
}
results.append(result)
return results
def _parse_md_tag(self, md_tag: str) -> list[int]:
"""Parse MD tag to extract mismatch positions."""
positions = []
read_pos = 1
i = 0
while i < len(md_tag):
if md_tag[i].isdigit():
num_str = ""
while i < len(md_tag) and md_tag[i].isdigit():
num_str += md_tag[i]
i += 1
if num_str:
read_pos += int(num_str)
elif md_tag[i] == "^":
i += 1
while i < len(md_tag) and md_tag[i].isalpha():
i += 1
elif md_tag[i].isalpha():
positions.append(read_pos)
read_pos += 1
i += 1
else:
i += 1
return positions
def _calculate_offtarget_score(self, mismatch_positions: list[int]) -> float:
"""Calculate off-target score based on mismatch count and positions.
Scoring principles (based on siRNA literature):
1. Total mismatch count - most fundamental property
2. Seed region (positions 2-8) - critical for target recognition
3. Position-specific weights - 5' end more important than 3' end
4. Continuous mismatches - clusters reduce binding more than scattered
Returns:
float: Off-target penalty score (lower = more likely off-target effect)
0.0 = perfect match (highest risk)
Higher scores = more mismatches (lower risk)
References: TODO: validate
- Jackson et al. 2003 (seed region importance)
- Birmingham et al. 2006 (position-specific effects)
- Huesken et al. 2005 (thermodynamic contributions)
"""
return _calculate_seed_offtarget_score(
mismatch_positions,
seed_start=self.seed_start,
seed_end=self.seed_end,
)
def _filter_and_rank(self, results: list[dict[str, Any]]) -> list[dict[str, Any]]:
"""Filter and rank results by off-target score."""
# as_score is None for records with no AS tag, so the dict default never fires.
results.sort(key=lambda x: (x["offtarget_score"], -int(x.get("as_score") or 0)))
return results
[docs]
class OffTargetAnalysisManager:
"""Manager class for comprehensive off-target analysis using BWA-MEM2."""
[docs]
def __init__(
self,
species: str,
transcriptome_path: str | Path | None = None,
mirna_path: str | Path | None = None,
transcriptome_index: str | Path | None = None,
mirna_index: str | Path | None = None,
):
"""Initialize the off-target analysis manager."""
self.species = species
self.transcriptome_path = Path(transcriptome_path) if transcriptome_path is not None else None
self.mirna_path = Path(mirna_path) if mirna_path is not None else None
self.transcriptome_index = Path(transcriptome_index) if transcriptome_index is not None else None
self.mirna_index = Path(mirna_index) if mirna_index is not None else None
[docs]
def analyze_mirna_off_targets(
self,
sequences: dict[str, str] | str | Path,
output_prefix: str | Path,
) -> tuple[Path, Path]:
"""Analyze miRNA off-targets using BWA-MEM2 in miRNA seed mode."""
if not self.mirna_index:
raise ValueError("miRNA index not provided")
if isinstance(sequences, str | Path):
sequences = parse_fasta_file(sequences)
analyzer = BwaAnalyzer(self.mirna_index, mode="mirna_seed")
results = analyzer.analyze_sequences(sequences)
output_path = Path(output_prefix)
tsv_path = output_path.parent / f"{output_path.name}_mirna_hits.tsv"
json_path = output_path.parent / f"{output_path.name}_mirna_hits.json"
self._write_mirna_results(results, tsv_path, json_path)
return tsv_path, json_path
[docs]
def analyze_transcriptome_off_targets(
self,
sequences: dict[str, str] | str | Path,
output_prefix: str | Path,
) -> tuple[Path, Path]:
"""Analyze transcriptome off-targets using BWA-MEM2 in transcriptome mode."""
if not self.transcriptome_index:
raise ValueError("Transcriptome index not provided")
if isinstance(sequences, str | Path):
sequences = parse_fasta_file(sequences)
analyzer = BwaAnalyzer(self.transcriptome_index, mode="transcriptome")
results = analyzer.analyze_sequences(sequences)
output_path = Path(output_prefix)
tsv_path = output_path.parent / f"{output_path.name}_transcriptome_hits.tsv"
json_path = output_path.parent / f"{output_path.name}_transcriptome_hits.json"
self._write_transcriptome_results(results, tsv_path, json_path)
return tsv_path, json_path
[docs]
def analyze_sirna_candidate(self, candidate: SiRNACandidate) -> dict[str, Any]:
"""Analyze a single siRNA candidate for off-targets."""
sequences = {candidate.id: candidate.guide_sequence}
results: dict[str, Any] = {
"candidate_id": candidate.id,
"guide_sequence": candidate.guide_sequence,
"mirna_hits": [],
"transcriptome_hits": [],
}
if self.mirna_index:
mirna_analyzer = BwaAnalyzer(self.mirna_index, mode="mirna_seed")
results["mirna_hits"] = mirna_analyzer.analyze_sequences(sequences)
if self.transcriptome_index:
transcriptome_analyzer = BwaAnalyzer(self.transcriptome_index, mode="transcriptome")
results["transcriptome_hits"] = transcriptome_analyzer.analyze_sequences(sequences)
return results
def _write_mirna_results(self, results: list[dict[str, Any]], tsv_path: str | Path, json_path: str | Path) -> None:
"""Write miRNA analysis results."""
# Write TSV
species_label = getattr(self, "species", "unknown")
enriched_results: list[dict[str, Any]] = []
with Path(tsv_path).open("w") as f:
f.write(
"qname\tqseq\tspecies\trname\tcoord\tstrand\tcigar\tmapq\tas_score\t"
"nm\tseed_mismatches\tofftarget_score\n"
)
for result in results:
row_species = result.get("species", species_label)
enriched = {
**result,
"species": row_species,
}
enriched_results.append(enriched)
f.write(
f"{enriched['qname']}\t{enriched['qseq']}\t{enriched['species']}\t{enriched['rname']}\t"
f"{enriched['coord']}\t{enriched['strand']}\t{enriched['cigar']}\t"
f"{enriched['mapq']}\t{enriched.get('as_score', 'NA')}\t{enriched['nm']}\t"
f"{enriched['seed_mismatches']}\t{enriched['offtarget_score']}\n"
)
# Write JSON
with Path(json_path).open("w") as f:
json.dump(enriched_results, f, indent=2)
def _write_transcriptome_results(
self, results: list[dict[str, Any]], tsv_path: str | Path, json_path: str | Path
) -> None:
"""Write transcriptome analysis results."""
# Write TSV
with Path(tsv_path).open("w") as f:
f.write("qname\tqseq\trname\tcoord\tstrand\tcigar\tmapq\tas_score\tnm\tseed_mismatches\tofftarget_score\n")
for result in results:
f.write(
f"{result['qname']}\t{result['qseq']}\t{result['rname']}\t"
f"{result['coord']}\t{result['strand']}\t{result['cigar']}\t"
f"{result['mapq']}\t{result.get('as_score', 'NA')}\t{result['nm']}\t"
f"{result['seed_mismatches']}\t{result['offtarget_score']}\n"
)
# Write JSON
with Path(json_path).open("w") as f:
json.dump(results, f, indent=2)
# =============================================================================
# Utility Functions
# =============================================================================
[docs]
def create_temp_fasta(sequences: dict[str, str]) -> str:
"""Create temporary FASTA file from sequences."""
with tempfile.NamedTemporaryFile(mode="w", suffix=".fasta", delete=False) as tmp_file:
temp_path = tmp_file.name
FastaUtils.write_dict_to_fasta(sequences, temp_path)
return temp_path
[docs]
def validate_and_write_sequences(
input_file: str, output_file: str, expected_length: int = 21
) -> tuple[int, int, list[str]]:
"""Validate siRNA sequences and write valid ones to output file."""
sequences = FastaUtils.parse_fasta_to_dict(input_file)
try:
valid_sequences = FastaUtils.validate_sirna_sequences(sequences, expected_length)
if valid_sequences:
FastaUtils.write_dict_to_fasta(valid_sequences, output_file)
else:
Path(output_file).touch()
invalid_count = len(sequences) - len(valid_sequences)
issues = [
f"{name}: Invalid (length={len(seq)}, expected={expected_length})"
for name, seq in sequences.items()
if name not in valid_sequences
]
return len(valid_sequences), invalid_count, issues
except ValueError as e:
Path(output_file).touch()
return 0, len(sequences), [str(e)]
[docs]
def build_bwa_index(fasta_file: str | Path, index_prefix: str | Path) -> Path:
"""Build BWA-MEM2 index for both transcriptome and miRNA off-target analysis."""
fasta_path = Path(fasta_file)
index_prefix_path = Path(index_prefix)
logger.info(f"Building BWA-MEM2 index from {fasta_path} with prefix {index_prefix_path}")
if not fasta_path.exists():
raise FileNotFoundError(f"Input FASTA file not found: {fasta_path}")
index_prefix_path.parent.mkdir(parents=True, exist_ok=True)
# Get absolute path to bwa-mem2 executable
bwa_path = _get_executable_path("bwa-mem2")
if not bwa_path:
raise FileNotFoundError("bwa-mem2 executable not found in PATH")
cmd = [bwa_path, "index", "-p", str(index_prefix_path), str(fasta_path)]
_validate_command_args(cmd)
try:
subprocess.run(cmd, capture_output=True, text=True, check=True, timeout=7200) # nosec B603
if not validate_index_files(index_prefix_path, "bwa-mem2"):
raise RuntimeError(
"BWA-MEM2 index files were not created correctly. This usually happens when the build "
"process runs out of memory (human transcriptomes can require 32GB+). Increase available RAM "
"or pre-build the index in an environment with more memory and retry."
)
logger.info(f"BWA-MEM2 index built successfully: {index_prefix_path}")
return index_prefix_path
except subprocess.CalledProcessError as e:
logger.error(f"BWA-MEM2 index build failed: {e.stderr}")
raise
except subprocess.TimeoutExpired:
logger.error("BWA-MEM2 index build timed out")
raise
[docs]
def validate_sirna_sequences(
sequences: dict[str, str], expected_length: int = 21
) -> tuple[dict[str, str], dict[str, str], list[str]]:
"""Validate siRNA sequences using existing FastaUtils."""
try:
valid_sequences = FastaUtils.validate_sirna_sequences(sequences, expected_length)
invalid_sequences = {name: seq for name, seq in sequences.items() if name not in valid_sequences}
issues = [
f"{name}: Invalid sequence (length={len(seq)}, expected={expected_length})"
for name, seq in invalid_sequences.items()
]
return valid_sequences, invalid_sequences, issues
except ValueError as e:
return {}, sequences, [str(e)]
[docs]
def parse_fasta_file(fasta_file: str | Path) -> dict[str, str]:
"""Parse FASTA file using existing FastaUtils."""
return FastaUtils.parse_fasta_to_dict(fasta_file)
[docs]
def write_fasta_file(sequences: dict[str, str], output_file: str) -> None:
"""Write sequences to FASTA file using existing FastaUtils."""
FastaUtils.write_dict_to_fasta(sequences, output_file)
[docs]
def validate_index_files(index_prefix: str | Path, tool: str = "bwa") -> bool:
"""Validate that index files exist for given tool."""
index_path = Path(index_prefix)
if tool in ("bwa", "bwa-mem2"):
required_extensions = [".amb", ".ann", ".bwt.2bit.64", ".pac"]
else:
logger.warning(f"Unknown tool for index validation: {tool}")
return False
for ext in required_extensions:
candidate = index_path.parent / f"{index_path.name}{ext}"
if not candidate.exists():
logger.debug(f"Missing index file: {candidate.name}")
return False
if candidate.stat().st_size == 0:
logger.debug(f"Index file is empty: {candidate}")
return False
return True
# =============================================================================
# Nextflow Integration Functions
# =============================================================================
def run_mirna_analysis_for_nextflow(
species: str,
sequences_file: str,
mirna_index: str | Path,
output_prefix: str | Path,
) -> tuple[str, str, str]:
"""Nextflow-compatible function for miRNA analysis."""
manager = OffTargetAnalysisManager(species=species, mirna_index=mirna_index)
output_root = Path(output_prefix)
try:
tsv_path, json_path = manager.analyze_mirna_off_targets(sequences_file, output_root)
# Create summary
summary_path = output_root.parent / f"{output_root.name}_mirna_summary.txt"
with summary_path.open("w") as f:
with json_path.open() as jf:
results = json.load(jf)
f.write(f"Species: {species}\n")
f.write(f"Total miRNA hits: {len(results)}\n")
f.write("Analysis completed successfully\n")
return str(tsv_path), str(json_path), str(summary_path)
except Exception as e:
error_summary = output_root.parent / f"{output_root.name}_mirna_error.txt"
with error_summary.open("w") as f:
f.write(f"miRNA analysis failed: {str(e)}\n")
return "", "", str(error_summary)
def run_transcriptome_analysis_for_nextflow(
species: str,
sequences_file: str,
transcriptome_index: str | Path,
output_prefix: str | Path,
) -> tuple[str, str, str]:
"""Nextflow-compatible function for transcriptome analysis."""
manager = OffTargetAnalysisManager(species=species, transcriptome_index=transcriptome_index)
output_root = Path(output_prefix)
try:
tsv_path, json_path = manager.analyze_transcriptome_off_targets(sequences_file, output_root)
# Create summary
summary_path = output_root.parent / f"{output_root.name}_transcriptome_summary.txt"
with summary_path.open("w") as f:
with json_path.open() as jf:
results = json.load(jf)
f.write(f"Species: {species}\n")
f.write(f"Total transcriptome hits: {len(results)}\n")
f.write("Analysis completed successfully\n")
return str(tsv_path), str(json_path), str(summary_path)
except Exception as e:
error_summary = output_root.parent / f"{output_root.name}_transcriptome_error.txt"
with error_summary.open("w") as f:
f.write(f"Transcriptome analysis failed: {str(e)}\n")
return "", "", str(error_summary)
def run_comprehensive_offtarget_analysis(
species: str,
sequences_file: str,
index_path: str,
output_prefix: str | Path,
mode: str = "transcriptome",
bwa_k: int = 12,
bwa_T: int = 15,
max_hits: int | None = None,
seed_start: int = 2,
seed_end: int = 8,
) -> tuple[str, str, str]:
"""Run comprehensive off-target analysis for Nextflow integration."""
output_root = Path(output_prefix)
try:
sequences = parse_fasta_file(sequences_file)
# Use BWA analyzer for comprehensive analysis
analyzer = BwaAnalyzer(
index_prefix=index_path,
mode=mode,
seed_length=bwa_k,
min_score=bwa_T,
max_hits=max_hits,
seed_start=seed_start,
seed_end=seed_end,
)
results = analyzer.analyze_sequences(sequences)
# Write results using pandas (much faster than manual loop)
tsv_path = output_root.parent / f"{output_root.name}.tsv"
json_path = output_root.parent / f"{output_root.name}.json"
summary_path = output_root.parent / f"{output_root.name}_summary.txt"
# Convert dict results to DataFrame
df = pd.DataFrame(results)
if df.empty:
# Create empty DataFrame with required schema columns (including species)
df = pd.DataFrame(columns=list(GenomeAlignmentSchema.__annotations__.keys()))
elif "species" not in df.columns:
df["species"] = species
else:
df["species"] = df["species"].fillna(species)
# Validate with Pandera schema
df = GenomeAlignmentSchema.validate(df, lazy=True)
# Write TSV and JSON (pandas handles efficiently)
df.to_csv(tsv_path, sep="\t", index=False)
df.to_json(json_path, orient="records", indent=2)
# Write summary
with summary_path.open("w") as f:
f.write(f"Species: {species}\n")
f.write(f"Total sequences analyzed: {len(sequences)}\n")
f.write(f"Total off-target hits: {len(results)}\n")
f.write(f"Analysis mode: {mode}\n")
f.write(f"Analysis parameters: bwa_k={bwa_k}, bwa_T={bwa_T}, max_hits={max_hits}\n")
f.write(f"Seed region: {seed_start}-{seed_end}\n")
f.write("Analysis completed successfully\n")
return str(tsv_path), str(json_path), str(summary_path)
except Exception as e:
error_summary = output_root.parent / f"{output_root.name}_error.txt"
with error_summary.open("w") as f:
f.write(f"Comprehensive off-target analysis failed: {str(e)}\n")
return "", "", str(error_summary)
[docs]
def run_bwa_alignment_analysis(
candidates_file: str | Path,
index_prefix: str | Path,
species: str,
output_dir: str | Path,
max_hits: int | None = None,
bwa_k: int = 12,
bwa_T: int = 15,
seed_start: int = 2,
seed_end: int = 8,
) -> Path:
"""Run BWA-MEM2 alignment analysis for candidate sequences using Pydantic models.
This is the main function called by OFFTARGET_ANALYSIS Nextflow module.
Args:
candidates_file: Path to FASTA file with candidate sequences
index_prefix: Path to BWA-MEM2 index prefix
species: Species identifier
output_dir: Directory to write results
max_hits: Maximum hits to report per candidate (``None`` = no limit / exhaustive)
bwa_k: BWA seed length parameter
bwa_T: BWA minimum score threshold
seed_start: Seed region start position (1-based)
seed_end: Seed region end position (1-based)
Returns:
Path to output directory containing results
"""
output_path = Path(output_dir)
output_path.mkdir(parents=True, exist_ok=True)
# Parse input sequences
sequences = parse_fasta_file(candidates_file)
# Determine candidate ID from filename (e.g., "candidate_0001.fasta" -> "candidate_0001")
candidate_id = Path(candidates_file).stem
# Create unique output prefix for this candidate-species combination
output_prefix = output_path / f"{candidate_id}_{species}"
# Run BWA-MEM2 analysis
analyzer = BwaAnalyzer(
index_prefix=index_prefix,
mode="transcriptome", # Always use transcriptome mode for genome analysis
seed_length=bwa_k,
min_score=bwa_T,
max_hits=max_hits,
seed_start=seed_start,
seed_end=seed_end,
)
results_dicts = analyzer.analyze_sequences(sequences)
# Convert dict results to OffTargetHit objects with validation
all_hits: list[OffTargetHit] = []
for hit_dict in results_dicts:
try:
# Parse coord string "chr1:12345" into integer
coord_str = hit_dict["coord"]
coord_int = int(coord_str.split(":")[1]) if ":" in coord_str else int(coord_str)
offtarget_hit = OffTargetHit(
qname=hit_dict["qname"],
qseq=hit_dict["qseq"],
species=species,
rname=hit_dict["rname"],
coord=coord_int,
strand=AlignmentStrand(hit_dict["strand"]),
cigar=hit_dict["cigar"],
mapq=hit_dict["mapq"],
as_score=hit_dict.get("as_score"),
nm=hit_dict["nm"],
seed_mismatches=hit_dict["seed_mismatches"],
offtarget_score=hit_dict["offtarget_score"],
)
all_hits.append(offtarget_hit)
except Exception as e:
logger.warning(f"Failed to validate off-target hit: {e}, skipping")
continue
# Write TSV analysis file using Pydantic models
analysis_file = Path(f"{output_prefix}_analysis.tsv")
with analysis_file.open("w") as f:
f.write(OffTargetHit.tsv_header() + "\n")
for hit in all_hits:
f.write(hit.to_tsv_row() + "\n")
# Write JSON file with validated data
json_file = Path(f"{output_prefix}_hits.json")
with json_file.open("w") as f:
json.dump([hit.model_dump() for hit in all_hits], f, indent=2)
# Derive summary statistics for alignment metrics
mean_mapq = statistics.fmean(hit.mapq for hit in all_hits) if all_hits else None
mean_mismatches = statistics.fmean(hit.nm for hit in all_hits) if all_hits else None
mean_seed_mismatches = statistics.fmean(hit.seed_mismatches for hit in all_hits) if all_hits else None
# Create validated summary using Pydantic model
summary = AnalysisSummary(
candidate_id=candidate_id,
species=species,
mode=AnalysisMode.TRANSCRIPTOME,
total_sequences=len(sequences),
total_hits=len(all_hits),
mean_mapq=mean_mapq,
mean_mismatches=mean_mismatches,
mean_seed_mismatches=mean_seed_mismatches,
)
# Write summary JSON file
summary_file = Path(f"{output_prefix}_summary.json")
with summary_file.open("w") as f:
# Add parameters to the output
summary_dict = summary.model_dump()
summary_dict["parameters"] = {
"bwa_k": bwa_k,
"bwa_T": bwa_T,
"max_hits": max_hits,
"seed_start": seed_start,
"seed_end": seed_end,
}
json.dump(summary_dict, f, indent=2)
logger.info(f"BWA analysis completed for {candidate_id} vs {species}: {len(all_hits)} hits")
return output_path
[docs]
def aggregate_offtarget_results( # noqa: PLR0912
results_dir: str | Path,
output_dir: str | Path,
genome_species: str,
) -> Path:
"""Aggregate transcriptome off-target analysis results using Pandera.
Uses pandas + Pandera for efficient bulk reading and validation instead of
manual line-by-line parsing with Pydantic models.
NOTE: This function ONLY aggregates genome/transcriptome hits. miRNA results
are aggregated separately by aggregate_mirna_results() to keep output files
distinct and properly typed.
Args:
results_dir: Directory containing individual analysis results
output_dir: Directory to write aggregated results
genome_species: Comma-separated list of genome species analyzed
Returns:
Path to output directory containing aggregated results
"""
results_path = Path(results_dir)
output_path = Path(output_dir)
output_path.mkdir(parents=True, exist_ok=True)
species_list = [s.strip() for s in genome_species.split(",") if s.strip()]
species_file_counts: dict[str, int] = {}
missing_species: list[str] = []
for species in species_list:
species_dir = results_path / species
count = len(list(species_dir.glob("*_analysis.tsv"))) if species_dir.exists() else 0
species_file_counts[species] = count
if count == 0:
missing_species.append(species)
# Collect ONLY genome/transcriptome TSV analysis files
# miRNA files are handled separately by aggregate_mirna_results()
analysis_files = list(results_path.glob("**/*_analysis.tsv"))
# Filter out miRNA files explicitly to avoid schema validation errors
analysis_files = [f for f in analysis_files if "mirna" not in f.name.lower()]
logger.info(f"Found {len(analysis_files)} transcriptome analysis files to aggregate")
if analysis_files:
# Read all files into DataFrames and concatenate (vectorized operation)
dfs = []
for analysis_file in analysis_files:
try:
# Pandas reads TSV much faster than manual line splitting
df = pd.read_csv(analysis_file, sep="\t")
# Validate schema with Pandera
df = GenomeAlignmentSchema.validate(df, lazy=True)
dfs.append(df)
except Exception as e:
logger.warning(f"Failed to read/validate {analysis_file}: {e}")
continue
# Concatenate all DataFrames at once (much faster than append in loop)
if dfs:
combined_df = pd.concat(dfs, ignore_index=True)
else:
# Create empty DataFrame with correct schema
combined_df = pd.DataFrame(columns=list(GenomeAlignmentSchema.__annotations__.keys()))
else:
# No files found - create empty DataFrame
combined_df = pd.DataFrame(columns=list(GenomeAlignmentSchema.__annotations__.keys()))
# Write combined results (pandas is much faster than manual TSV writing)
combined_tsv = output_path / "combined_offtargets.tsv"
combined_df.to_csv(combined_tsv, sep="\t", index=False)
# Write JSON (pandas handles serialization)
combined_json = output_path / "combined_offtargets.json"
combined_df.to_json(combined_json, orient="records", indent=2)
species_counts = _compute_species_counts(combined_df)
if not species_counts:
species_counts = dict.fromkeys(species_list, 0) if species_list else {}
for species in species_list:
species_counts.setdefault(species, 0)
human_hits, other_hits = human_vs_other_totals(species_counts)
logger.info(
f"Aggregated {len(combined_df)} transcriptome off-target hits from {len(analysis_files)} files using pandas"
)
# Prepare summary metadata
summary_json = output_path / "combined_summary.json"
# Create validated aggregated summary
summary_status = "completed" if not missing_species else "partial"
summary = AggregatedOffTargetSummary(
species_analyzed=species_list,
analysis_files_processed=len(analysis_files),
total_results=len(combined_df),
combined_tsv=combined_tsv,
combined_json=combined_json,
summary_file=summary_json,
hits_per_species=species_counts,
human_hits=human_hits,
other_species_hits=other_hits,
species_file_counts=species_file_counts,
missing_species=missing_species,
status=summary_status,
)
# Write summary JSON
with summary_json.open("w") as f:
json.dump(summary.model_dump(mode="json"), f, indent=2)
# Write final summary text file
final_summary = output_path / "final_summary.txt"
with final_summary.open("w") as f:
f.write("Off-Target Analysis Aggregation Summary\n")
f.write("=" * 50 + "\n\n")
# Check if transcriptome analysis was performed
if len(analysis_files) == 0:
f.write("TRANSCRIPTOME ANALYSIS STATUS: NOT PERFORMED\n")
f.write("-" * 50 + "\n")
f.write("Reason: No transcriptome FASTAs or BWA indices were provided.\n")
f.write("Result: Only lightweight miRNA seed match analysis was run.\n\n")
f.write("To enable transcriptome off-target analysis:\n")
f.write(" • Provide --genome_fastas (transcriptome) 'species:path,species2:path2'\n")
f.write(" OR\n")
f.write(" • Provide --genome_indices 'species:index,species2:index2'\n\n")
f.write("=" * 50 + "\n\n")
if missing_species:
warning_list = ", ".join(missing_species)
f.write("WARNINGS\n")
f.write("-" * 50 + "\n")
f.write(
"No transcriptome alignment files were produced for the following species: "
f"{warning_list}. This usually indicates the BWA-MEM2 indexing stage ran out of memory.\n"
)
f.write(
"Increase --max_memory (32GB+ recommended for human transcriptomes) or pre-build indices on a host with more RAM.\n\n"
)
# Results summary
f.write("RESULTS SUMMARY\n")
f.write("-" * 50 + "\n")
f.write(f"Transcriptome off-target hits: {len(combined_df)}\n")
f.write(f"Human hits: {human_hits}\n")
f.write(f"Other species hits: {other_hits}\n")
# Show species list or note if empty
if species_list:
f.write(f"Species requested for analysis: {', '.join(species_list)}\n")
else:
f.write("Species requested for analysis: (none - miRNA-only mode)\n")
if species_counts:
f.write("Per-species hit counts:\n")
for species, count in sorted(species_counts.items()):
f.write(f" {species}: {count}\n")
f.write(f"Transcriptome analysis files processed: {len(analysis_files)}\n\n")
# Explain what the output files contain
f.write("OUTPUT FILES\n")
f.write("-" * 50 + "\n")
if len(combined_df) == 0:
f.write(f"• {combined_tsv.name}: Header only (no hits found)\n")
f.write(f"• {combined_json.name}: Empty array (no hits found)\n")
f.write(f"• {summary_json.name}: Metadata only\n\n")
f.write(
"Note: Empty data files indicate NO problematic transcriptome off-targets were detected - this is GOOD!\n"
)
f.write("Your siRNA candidates are clean at the transcriptome alignment level.\n\n")
f.write("For miRNA seed match analysis results, see:\n")
f.write(" ../mirna/mirna_analysis.tsv\n")
f.write(" ../mirna/mirna_summary.json\n")
else:
f.write(f"• {combined_tsv.name}: {len(combined_df)} off-target hits (TSV format)\n")
f.write(f"• {combined_json.name}: {len(combined_df)} off-target hits (JSON format)\n")
f.write(f"• {summary_json.name}: Analysis metadata and statistics\n")
if missing_species:
logger.warning(
"Transcriptome aggregation completed with missing species: %s. "
"Likely cause: insufficient memory while building BWA-MEM2 indices.",
", ".join(missing_species),
)
else:
logger.info(f"Wrote aggregated results to {output_path}")
return output_path
[docs]
def run_mirna_seed_analysis(
candidates_file: str | Path,
candidate_id: str,
mirna_db: str, # Review, can this be linked to a class describing all miRNA database protocol/ABC?
mirna_species: list[str],
output_dir: str | Path,
backend: MiRNASeedBackend | str = MiRNASeedBackend.PYAHOCORASICK,
seed_start: int = 2,
seed_end: int = 8,
) -> Path:
"""Run miRNA seed match analysis for candidate sequences.
This function uses the MiRNADatabaseManager to download and cache miRNA databases,
builds BWA indices if needed, and performs seed match analysis.
The scan produces *raw* alignments: the guide seed window placed at every position
along each miRNA. Only alignments where the guide seed lands on the miRNA's own seed
region (0-based ``coord == seed_start - 1``) are counted as *hits* in the filtered
outputs and summary ``total_hits``; perfect matches in non-seed regions are retained
in the ``*_raw`` files but are not real miRNA seed off-targets.
Args:
candidates_file: Path to FASTA file with candidate sequences
candidate_id: Candidate identifier
mirna_db: miRNA database name (mirgenedb, mirbase, etc.)
mirna_species: List of species to analyze against
output_dir: Directory to write results
backend: miRNA seed backend to use for analysis (pyahocorasick by default)
seed_start: Seed region start position (1-based, default 2)
seed_end: Seed region end position (1-based, default 8)
Returns:
Path to output directory containing results
"""
output_path = Path(output_dir)
output_path.mkdir(parents=True, exist_ok=True)
# Initialize miRNA database manager
manager = MiRNADatabaseManager()
# Parse input sequences
sequences = parse_fasta_file(candidates_file)
all_raw_hits = [] # All raw alignments from BWA
species_raw_stats = {}
species_filtered_stats = {}
resolved_backend = MiRNASeedBackend(backend)
logger.info(f"Running miRNA seed analysis for {candidate_id}")
logger.info(f"Database: {mirna_db}, Species: {mirna_species}, Backend: {resolved_backend.value}")
for species in mirna_species:
try:
# Get or download miRNA database for this species
logger.info(f"Processing miRNA database for species: {species}")
db_fasta_path = manager.get_database(mirna_db, species)
if db_fasta_path is None or not db_fasta_path.exists():
logger.warning(f"miRNA database not available for {species}, skipping")
continue
if resolved_backend == MiRNASeedBackend.BWA:
mirna_sequences = parse_fasta_file(db_fasta_path)
index_prefix = db_fasta_path.with_suffix("")
if not validate_index_files(index_prefix, "bwa-mem2"):
logger.info(f"Building BWA index for {species} miRNA database")
build_bwa_index(db_fasta_path, index_prefix)
analyzer = BwaAnalyzer(
index_prefix=index_prefix,
mode="mirna_seed",
seed_length=6,
min_score=6,
max_hits=_mirna_max_hits(),
seed_start=seed_start,
seed_end=seed_end,
)
bwa_results = analyzer.analyze_sequences(sequences)
results = _normalize_bwa_mirna_seed_hits(
bwa_results,
sequences=sequences,
mirna_sequences=mirna_sequences,
seed_start=seed_start,
seed_end=seed_end,
max_mismatches=2,
)
else:
mirna_sequences = parse_fasta_file(db_fasta_path)
results = scan_mirna_seed_matches(
sequences,
mirna_sequences,
backend=resolved_backend,
seed_start=seed_start,
seed_end=seed_end,
max_mismatches=2,
max_hits=_mirna_max_hits(),
)
results_df = _build_mirna_alignment_frame(
results,
species=species,
database=mirna_db,
)
# Validate and coerce types using Pandera schema
try:
validated_df = MiRNAAlignmentSchema.validate(results_df, lazy=True)
all_raw_hits.append(validated_df)
species_raw_stats[species] = len(validated_df)
logger.info(f"Species {species}: {len(validated_df)} miRNA alignments validated")
except Exception as validation_error:
logger.error(f"Failed to validate miRNA hits for {species}: {validation_error}")
species_raw_stats[species] = 0
species_filtered_stats[species] = 0
continue
except _MiRNASeedBackendUnavailableError as backend_error:
message = (
f"miRNA seed analysis backend '{resolved_backend.value}' is unavailable for species "
f"{species}: {backend_error}"
)
logger.error(message)
raise RuntimeError(message) from backend_error
except Exception as e:
logger.error(f"Failed to process miRNA analysis for {species}: {e}")
species_raw_stats[species] = 0
species_filtered_stats[species] = 0
# Concatenate all validated DataFrames from different species
if all_raw_hits:
df_raw = pd.concat(all_raw_hits, ignore_index=True)
# A raw alignment is a real miRNA seed off-target only when the guide seed lands
# on the miRNA's own seed region, i.e. the 0-based alignment start equals the seed
# start offset (seed_start - 1). Perfect matches elsewhere in the miRNA (e.g. its 3'
# region) are retained in the *_raw outputs but are NOT counted as seed hits.
seed_region_coord = seed_start - 1
df_filtered = df_raw[df_raw["coord"] == seed_region_coord].reset_index(drop=True)
# Calculate per-species filtered stats
for species in mirna_species:
species_filtered_stats[species] = len(df_filtered[df_filtered["species"] == species])
# Write RAW hits TSV (all alignments)
raw_analysis_file = output_path / f"{candidate_id}_mirna_analysis_raw.tsv"
df_raw.to_csv(raw_analysis_file, sep="\t", index=False)
# Write FILTERED hits TSV (quality-filtered)
filtered_analysis_file = output_path / f"{candidate_id}_mirna_analysis.tsv"
df_filtered.to_csv(filtered_analysis_file, sep="\t", index=False)
# Write raw hits JSON
raw_json_file = output_path / f"{candidate_id}_mirna_hits_raw.json"
df_raw.to_json(raw_json_file, orient="records", indent=2)
# Write filtered hits JSON
filtered_json_file = output_path / f"{candidate_id}_mirna_hits.json"
df_filtered.to_json(filtered_json_file, orient="records", indent=2)
total_filtered = len(df_filtered)
total_raw = len(df_raw)
else:
# No hits - create empty DataFrame with proper schema columns
df_empty = pd.DataFrame(columns=list(MiRNAAlignmentSchema.to_schema().columns.keys()))
raw_analysis_file = output_path / f"{candidate_id}_mirna_analysis_raw.tsv"
df_empty.to_csv(raw_analysis_file, sep="\t", index=False)
filtered_analysis_file = output_path / f"{candidate_id}_mirna_analysis.tsv"
df_empty.to_csv(filtered_analysis_file, sep="\t", index=False)
raw_json_file = output_path / f"{candidate_id}_mirna_hits_raw.json"
df_empty.to_json(raw_json_file, orient="records", indent=2)
filtered_json_file = output_path / f"{candidate_id}_mirna_hits.json"
df_empty.to_json(filtered_json_file, orient="records", indent=2)
total_filtered = 0
total_raw = 0
# Create validated summary using Pydantic model
summary = MiRNASummary(
candidate_id=candidate_id,
mirna_database=mirna_db,
species_analyzed=mirna_species,
total_sequences=len(sequences),
total_hits=total_filtered, # Filtered hits count
hits_per_species=species_raw_stats, # Raw hits per species
total_raw_alignments=total_raw, # Total raw alignments
)
# Write summary JSON file
summary_file = output_path / f"{candidate_id}_mirna_summary.json"
with summary_file.open("w") as f:
json.dump(summary.model_dump(mode="json"), f, indent=2)
logger.info(
f"miRNA seed analysis completed for {candidate_id}: "
f"{total_filtered} filtered high-quality matches "
f"(from {total_raw} raw alignments)"
)
return output_path
[docs]
def aggregate_mirna_results(
results_dir: str | Path,
output_dir: str | Path,
mirna_db: str,
mirna_species: str,
) -> Path:
"""Aggregate miRNA seed analysis results from multiple candidates using pandas.
Uses pandas + Pandera for efficient bulk reading and validation instead of
manual line-by-line parsing with Pydantic models.
Args:
results_dir: Directory containing individual miRNA analysis results
output_dir: Directory to write aggregated results
mirna_db: miRNA database used for analysis
mirna_species: Comma-separated list of species analyzed
Returns:
Path to output directory containing aggregated results
"""
results_path = Path(results_dir)
output_path = Path(output_dir)
output_path.mkdir(parents=True, exist_ok=True)
species_list = [s.strip() for s in mirna_species.split(",") if s.strip()]
# Collect all miRNA analysis files using pandas (much faster than manual parsing)
analysis_files = list(results_path.glob("**/*_mirna_analysis.tsv"))
if analysis_files:
# Read all files into DataFrames and concatenate (vectorized operation)
dfs = []
candidate_stats = {}
for analysis_file in analysis_files:
try:
# Extract candidate ID from filename
candidate_id = analysis_file.stem.replace("_mirna_analysis", "")
# Pandas reads TSV much faster than manual line splitting
df = pd.read_csv(analysis_file, sep="\t")
# Validate schema with Pandera
df = MiRNAAlignmentSchema.validate(df, lazy=True)
# Track hits per candidate
candidate_stats[candidate_id] = len(df)
dfs.append(df)
except Exception as e:
logger.warning(f"Failed to read/validate {analysis_file}: {e}")
continue
# Concatenate all DataFrames at once (much faster than append in loop)
if dfs:
combined_df = pd.concat(dfs, ignore_index=True)
else:
# Create empty DataFrame with correct schema
combined_df = pd.DataFrame(columns=list(MiRNAHit.model_fields.keys()))
else:
# No files found - create empty DataFrame
combined_df = pd.DataFrame(columns=list(MiRNAHit.model_fields.keys()))
candidate_stats = {}
# Write combined results (pandas is much faster than manual TSV writing)
combined_tsv = output_path / "combined_mirna_hits.tsv"
combined_df.to_csv(combined_tsv, sep="\t", index=False)
# Write JSON (pandas handles serialization)
combined_json = output_path / "combined_mirna_hits.json"
combined_df.to_json(combined_json, orient="records", indent=2)
logger.info(f"Aggregated {len(combined_df)} miRNA hits from {len(analysis_files)} files using pandas")
# Calculate statistics using pandas groupby (much faster than loops)
species_stats = _compute_species_counts(combined_df)
if not species_stats:
species_stats = dict.fromkeys(species_list, 0) if species_list else {}
for species in species_list:
species_stats.setdefault(species, 0)
human_hits, other_hits = human_vs_other_totals(species_stats)
# Create validated summary using Pydantic model
summary = AggregatedMiRNASummary(
total_mirna_hits=len(combined_df),
mirna_database=mirna_db,
species_analyzed=species_list,
hits_per_species=species_stats,
hits_per_candidate=candidate_stats,
analysis_files_processed=len(analysis_files),
total_candidates=len(candidate_stats),
combined_tsv=combined_tsv,
combined_json=combined_json,
summary_file=output_path / "combined_mirna_summary.json",
human_hits=human_hits,
other_species_hits=other_hits,
)
# Write summary JSON
summary_json = output_path / "combined_mirna_summary.json"
with summary_json.open("w") as f:
json.dump(summary.model_dump(mode="json"), f, indent=2)
# Write final summary text file
final_summary = output_path / "final_mirna_summary.txt"
with final_summary.open("w") as f:
f.write("miRNA Seed Match Analysis Aggregation Summary\n")
f.write("=" * 50 + "\n")
f.write(f"Total miRNA seed matches: {len(combined_df)}\n")
f.write(f"Database: {mirna_db}\n")
f.write(f"Species analyzed: {', '.join(species_list)}\n")
f.write(f"Candidates analyzed: {len(candidate_stats)}\n")
f.write(f"Analysis files processed: {len(analysis_files)}\n")
f.write(f"Human hits: {human_hits}\n")
f.write(f"Other species hits: {other_hits}\n")
f.write("\nHits per species:\n")
for species, count in species_stats.items():
f.write(f" {species}: {count}\n")
f.write("\nOutput files:\n")
f.write(f" - Combined TSV: {combined_tsv.name}\n")
f.write(f" - Combined JSON: {combined_json.name}\n")
f.write(f" - Summary JSON: {summary_json.name}\n")
logger.info(f"Wrote aggregated miRNA results to {output_path}")
return output_path
# Export all main functions and classes
__all__ = [
# Core classes
"BwaAnalyzer",
"OffTargetAnalysisManager",
"MiRNASeedBackend",
# Utility functions
"create_temp_fasta",
"validate_and_write_sequences",
"build_bwa_index",
"validate_sirna_sequences",
"parse_fasta_file",
"write_fasta_file",
"check_tool_availability",
"validate_index_files",
"mirna_seed_hit_identity",
"normalize_mirna_seed_hit",
# Nextflow integration functions (called directly from Nextflow modules)
"run_bwa_alignment_analysis",
"aggregate_offtarget_results",
"scan_mirna_seed_matches",
"run_mirna_seed_analysis",
"aggregate_mirna_results",
]