BioSight/python/phi/output.py

"""
phi.output — FASTQ, Adleman graph, and PCR primer protocol

Transforms Φ encoding records (from phi.embed.encode_phi) into
standard bioinformatics formats consumable by external tools.

Formats:
  FASTQ          — pipe through Jellyfish/KMC/BLAST
  Adleman graph  — feed into BioComputing's connectNodes()
  PCR protocol   — wet-lab filtering with allele-specific primers

Dependencies: phi.embed (for encode_phi output dict format).
"""

from __future__ import annotations

import math
from typing import Dict, List

# ── Primer sequences ─────────────────────────────────────────────────────

PRIMER_DESIGN = {
    "pas_primer": {
        "sequence": "CCCCCC",
        "target": "consistency_dna = GGGGGG (all 6 rules pass)",
        "tm": 36.0,
        "purpose": "Capture PASS strands — amplify consistent equations only",
    },
    "fail_primer": {
        "sequence": "AAAAAA",
        "target": "consistency_dna = TTTTTT (all 6 rules fail)",
        "tm": 12.0,
        "purpose": "Capture FAIL strands — quarantine inconsistent equations",
    },
    "universal_primer": {
        "sequence": "ABCABC",
        "target": "F_dna region (first 8 bases)",
        "tm": 15.0,
        "purpose": "Amplify all encoded equations regardless of consistency",
    },
}


# ── FASTQ ────────────────────────────────────────────────────────────────

def to_fastq(record: Dict) -> str:
    """Format a Φ encoding record as a 4-line FASTQ entry.

    FASTQ is the standard DNA sequencing format.  This lets us pipe
    equation encodings through any existing bioinformatics pipeline
    without tool modification.
    """
    eq = record["equation"]
    seq = record["dna_sequence"]
    qual = record["quality_scores"]
    n_cons = sum(1 for v in record["consistency"].values() if v)
    eq_hash = record["sha256_prefix"]
    pass_flag = "PASS" if record["consistency_pass"] else "FAIL"

    header = f"@{eq_hash} n_consistency={n_cons}/6 pass={pass_flag} len={len(eq)}"
    qual_line = qual * (len(seq) // len(qual)) if len(seq) > len(qual) else qual

    return f"{header}\n{seq}\n+\n{qual_line}\n"


# ── Adleman integration ──────────────────────────────────────────────────

def to_adleman_graph(records: List[Dict]) -> Dict:
    """Build an Adleman/Lipton-compatible graph from Φ-encoded equations.

    Each equation is a vertex.  Edges represent Fisher-distance similarity
    (d_F < 0.5) between their F(E) byte-class histograms.

    Feed the output into BioComputing's `hampath.connectNodes()` or
    `sattv.createNodes()` to generate wet-lab DNA sequences.
    """
    vertices = {}
    for i, rec in enumerate(records):
        vid = f"E{i}"
        vertices[vid] = {
            "equation": rec["equation"],
            "dna": rec["dna_sequence"],
            "consistency": rec["consistency_pass"],
        }

    edges = []
    vert_ids = list(vertices.keys())
    for i in range(len(vert_ids)):
        for j in range(i + 1, len(vert_ids)):
            Fi = records[i]["F"]
            Fj = records[j]["F"]
            dist = _fisher_distance(Fi, Fj)
            if dist < 0.5:
                edges.append((vert_ids[i], vert_ids[j], round(dist, 4)))

    return {
        "vertices": vertices,
        "edges": edges,
        "format": "adleman_graph_v1",
        "n_vertices": len(vertices),
        "n_edges": len(edges),
    }


# ── PCR protocol ─────────────────────────────────────────────────────────

def design_filtering_protocol(records: List[Dict]) -> Dict:
    """Design a PCR filtering protocol for a batch of encoded equations.

    Uses allele-specific PCR (Newton et al. 1989) — the 24°C Tm gap
    between PASS and FAIL primers provides single-base discrimination.
    """
    n_pass = sum(1 for r in records if r["consistency_pass"])
    n_fail = len(records) - n_pass

    return {
        "protocol": "allele_specific_pcr_v1",
        "reference": "Newton et al. (1989) Nucleic Acids Res 17:2503; "
                      "allele-specific PCR for single-base discrimination",
        "pas_primer": PRIMER_DESIGN["pas_primer"]["sequence"],
        "fail_primer": PRIMER_DESIGN["fail_primer"]["sequence"],
        "universal_primer": PRIMER_DESIGN["universal_primer"]["sequence"],
        "annealing_tm_pass": PRIMER_DESIGN["pas_primer"]["tm"],
        "annealing_tm_fail": PRIMER_DESIGN["fail_primer"]["tm"],
        "annealing_tm_universal": PRIMER_DESIGN["universal_primer"]["tm"],
        "expected_amplification_pass": n_pass,
        "expected_amplification_fail": n_fail,
    }


# ── Fisher distance ↓ (internal; no circular dep) ────────────────────────

def _fisher_distance(p: List[float], q: List[float]) -> float:
    """Fisher-Rao distance between two points on Δ₇: d_F = 2·acos(Σ√(p_i·q_i))."""
    dot = sum(math.sqrt(pi * qi) for pi, qi in zip(p, q))
    dot = max(-1.0, min(1.0, dot))
    return 2.0 * math.acos(dot)