"""Check the Cu/Zn ligands of SOD1-like proteins below a PAINT key-residue (IKR) loss.

PAINT's IKR at PTHR10003 node PTN008300295 (Diptera) negates superoxide dismutase
activity for Drosophila Sodq and its mosquito ortholog. An IKR claims that residues the
activity depends on are lost. This script tests that claim: it reads the metal-ligand
positions of human SOD1 (P00441) from UniProt's BINDING features, aligns each query to
SOD1, and reports the residue each query has at every ligand position. Drosophila Sod1
(P61851), an active dismutase, is the positive control.

Sequences and features are fetched from the UniProt REST API. Writes
``ikr_metal_ligands.md`` beside this script.

Usage:
    uv run python projects/IRD_EVIDENCE/ikr_metal_ligands.py
"""

from __future__ import annotations

import re
import sys
from pathlib import Path

import requests
from Bio import Align
from Bio.Align import substitution_matrices

HERE = Path(__file__).parent
OUT = HERE / "ikr_metal_ligands.md"
REFERENCE = "P00441"
QUERIES = {
    "P61851": "Drosophila Sod1 (active SOD1, positive control)",
    "Q9VBK8": "Drosophila Sodq (below the IKR)",
    "A0A1S4GW26": "Anopheles ortholog (below the IKR)",
}


def fetch_entry(accession: str) -> str:
    """Return the UniProt flat-file text for one accession."""
    response = requests.get(f"https://rest.uniprot.org/uniprotkb/{accession}.txt", timeout=60)
    response.raise_for_status()
    return response.text


def sequence_of(entry: str) -> str:
    """Return the sequence from a UniProt flat file."""
    block = entry.split("\nSQ   ", 1)[1].split("\n", 1)[1].split("//", 1)[0]
    return re.sub(r"\s", "", block)


def ligand_positions(entry: str) -> list[tuple[int, str]]:
    """Return (1-based position, ligand) for each metal BINDING feature.

    >>> ligand_positions('FT   BINDING         47\\nFT                   /ligand="Cu cation"\\n')
    [(47, 'Cu cation')]
    """
    found = re.findall(r'^FT   BINDING\s+(\d+)\n(?:FT {19}/.*\n)*?FT {19}/ligand="([^"]+)"', entry, re.M)
    return [(int(pos), ligand) for pos, ligand in found]


def aligned_residue(reference: str, query: str, position: int) -> str:
    """Residue the query aligns to the reference's 1-based ``position`` ('-' for a gap)."""
    aligner = Align.PairwiseAligner(mode="global")
    aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
    aligner.open_gap_score, aligner.extend_gap_score = -10, -0.5
    alignment = aligner.align(reference, query)[0]
    for (r_start, r_end), (q_start, _) in zip(*alignment.aligned):
        if r_start <= position - 1 < r_end:
            return query[q_start + position - 1 - r_start]
    return "-"


def main() -> int:
    """Write the ligand table."""
    ref_entry = fetch_entry(REFERENCE)
    reference = sequence_of(ref_entry)
    ligands = ligand_positions(ref_entry)
    header = "| Protein | " + " | ".join(
        f"{reference[p - 1]}{p} ({'Cu' if 'Cu' in lig else 'Zn'})" for p, lig in ligands
    ) + " | Ligands kept |"
    lines = [
        "# Cu/Zn ligands below the PTHR10003 PTN008300295 IKR",
        "",
        f"Ligand positions are human SOD1 ({REFERENCE}) UniProt BINDING features; each query is",
        "globally aligned to SOD1 (BLOSUM62, gap open -10, extend -0.5). Generated by",
        "`ikr_metal_ligands.py`.",
        "",
        header,
        "|" + "---|" * (len(ligands) + 2),
    ]
    for accession, label in QUERIES.items():
        query = sequence_of(fetch_entry(accession))
        residues = [aligned_residue(reference, query, p) for p, _ in ligands]
        kept = sum(r == reference[p - 1] for r, (p, _) in zip(residues, ligands))
        lines.append(f"| {accession} {label} | " + " | ".join(residues) + f" | {kept}/{len(ligands)} |")
    OUT.write_text("\n".join(lines) + "\n")
    print("\n".join(lines))
    return 0


if __name__ == "__main__":
    sys.exit(main())