SimurgArray/tests/oracle/test_p1_real_files.py

"""Phase P1 acceptance on real files (plan §18 P1): parsers agree field-by-field with Illumina's BeadArrayFiles
(BSD-2) and GTCs round-trip byte-identically. Needs ``scripts/fetch_data.sh p0``; GTC oracles are produced with
gtc2vcf ``idat2gtc`` (see docs/validation/P1.md).
"""

from __future__ import annotations

import numpy as np
import pytest
from tests.conftest import (
    BIOC_IDATS,
    GDA_EGT,
    GDA_SHEET,
    GSA_BPM,
    GSA_CSV,
    GSA_EGT,
    GSA_IDATS,
    ORACLE_GTCS,
    need,
)

from simurg.core.errors import FormatError
from simurg.formats.bpm import read_bpm
from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest
from simurg.formats.discovery import discover_idats
from simurg.formats.egt import read_egt
from simurg.formats.gtc import encode_gtc, read_gtc
from simurg.formats.idat import check_pair, read_idat
from simurg.formats.samplesheet import parse_samplesheet

pytestmark = pytest.mark.data


@pytest.fixture(scope="module")
def gsa_bpm():  # type: ignore[no-untyped-def]
    need(GSA_BPM)
    return read_bpm(GSA_BPM)


def test_gsa_idats_all_parse_and_pair() -> None:
    need(GSA_IDATS)
    d = discover_idats(GSA_IDATS)
    assert len(d.pairs) == 36 and all(p.complete for p in d.pairs.values())
    ref = None
    for pair in d.pairs.values():
        g, r = read_idat(pair.grn), read_idat(pair.red)  # type: ignore[arg-type]
        check_pair(g, r)
        assert g.barcode == pair.barcode and g.position == pair.position
        assert g.scan_time_iso  # every demo IDAT has a scan entry
        ids = np.sort(g.address_ids)
        if ref is None:
            ref = ids
        assert np.array_equal(ids, ref)


def test_bioc_idats() -> None:
    need(BIOC_IDATS)
    ok = {
        "4019585376_B_Red.idat": 381079,
        "4019585376_B_Red.idat.gz": 381079,
        "5723646052_R02C02_Grn.idat": 622399,
        "5723646052_R02C02_Grn.idat.gz": 622399,
    }
    for name, n in ok.items():
        assert read_idat(BIOC_IDATS / name).n_addresses == n
    for name in ("4343238080_A_Grn.idat", "6016741005_A_Grn.idat"):
        with pytest.raises(FormatError) as e:
            read_idat(BIOC_IDATS / name)
        assert e.value.code == "E_IDAT_ENCRYPTED"


def test_idat_addresses_cover_manifest(gsa_bpm) -> None:  # type: ignore[no-untyped-def]  # precursor of IN-027
    need(GSA_IDATS)
    pair = next(iter(discover_idats(GSA_IDATS).pairs.values()))
    ids = read_idat(pair.grn).address_ids  # type: ignore[arg-type]
    wanted = np.concatenate([gsa_bpm.address_a, gsa_bpm.address_b[gsa_bpm.address_b != 0]])
    assert np.isin(wanted, ids).mean() == 1.0


def test_gtc_matches_beadarrayfiles_and_roundtrips(beadarrayfiles) -> None:  # type: ignore[no-untyped-def]
    need(ORACLE_GTCS)
    files = sorted(ORACLE_GTCS.glob("*.gtc"))
    assert files
    GenotypeCalls = beadarrayfiles.GenotypeCalls
    for f in files:
        g = read_gtc(f)
        assert encode_gtc(g) == f.read_bytes(), f"round-trip differs: {f}"
        b = GenotypeCalls(str(f))
        assert np.array_equal(g.raw_x, b.get_raw_x_intensities())
        assert np.array_equal(g.raw_y, b.get_raw_y_intensities())
        assert g.genotypes.tolist() == list(b.get_genotypes())
        assert np.array_equal(g.genotype_scores, b.get_genotype_scores(), equal_nan=True)
        assert np.array_equal(g.b_allele_freqs, b.get_ballele_freqs(), equal_nan=True)
        assert np.array_equal(g.logr_ratios, b.get_logr_ratios(), equal_nan=True)
        assert g.base_call_strings() == [x.decode() for x in b.get_base_calls()]
        assert (g.call_rate, g.logr_dev, g.gc10, g.gc50) == (
            b.get_call_rate(),
            b.get_logr_dev(),
            b.get_gc10(),
            b.get_gc50(),
        )
        assert g.gender == b.get_gender()
        assert (g.num_calls, g.num_no_calls, g.num_intensity_only) == (
            b.get_num_calls(),
            b.get_num_no_calls(),
            b.get_num_intensity_only(),
        )
        assert g.percentiles_x.tolist() == list(b.get_percentiles_x())
        xf = [
            (t.offset_x, t.offset_y, t.scale_x, t.scale_y, t.shear, t.theta) for t in b.get_normalization_transforms()
        ]
        mine = [
            tuple(float(x[k]) for k in ("offset_x", "offset_y", "scale_x", "scale_y", "shear", "theta"))
            for x in g.normalization_transforms
        ]
        assert mine == [tuple(float(v) for v in t) for t in xf]


def test_gtc_norm_transform_count_matches_bpm_bins(gsa_bpm) -> None:  # type: ignore[no-untyped-def]
    need(ORACLE_GTCS)
    g = read_gtc(sorted(ORACLE_GTCS.glob("*.gtc"))[0])
    assert len(g.normalization_transforms) == len(np.unique(gsa_bpm.norm_bin_ids()))
    assert g.num_snps == gsa_bpm.num_loci


@pytest.mark.slow
def test_bpm_matches_beadarrayfiles(gsa_bpm, beadarrayfiles) -> None:  # type: ignore[no-untyped-def]
    b = beadarrayfiles.BeadPoolManifest(str(GSA_BPM))
    m = gsa_bpm
    assert m.names == b.names
    assert m.snps == b.snps
    assert m.chroms == b.chroms
    assert m.map_infos.tolist() == b.map_infos
    assert m.address_a.tolist() == b.addresses
    assert m.assay_types.tolist() == b.assay_types
    assert m.normalization_lookups().tolist() == b.normalization_lookups
    rs = beadarrayfiles.RefStrand
    assert [rs.from_string(s) for s in m.ref_strands] == b.ref_strands
    assert m.manifest_name == b.manifest_name and m.control_config == b.control_config


@pytest.mark.slow
@pytest.mark.parametrize("egt_path", [GSA_EGT, GDA_EGT], ids=["gsa", "gda"])
def test_egt_matches_beadarrayfiles(egt_path, beadarrayfiles) -> None:  # type: ignore[no-untyped-def]
    need(egt_path)
    e = read_egt(egt_path)
    with open(egt_path, "rb") as fh:
        b = beadarrayfiles.ClusterFile.read_cluster_file(fh)
    assert e.num_records == len(b.records) if hasattr(b, "records") else True
    idx = e.name_index()
    for name in e.names[:: max(1, e.num_records // 5000)]:  # ~5000 evenly spaced records
        rec = b.get_record(name)
        i = idx[name]
        for k, stats in enumerate((rec.aa_cluster_stats, rec.ab_cluster_stats, rec.bb_cluster_stats)):
            assert e.n[i, k] == stats.N
            assert e.theta_mean[i, k] == np.float32(stats.theta_mean)
            assert e.theta_dev[i, k] == np.float32(stats.theta_dev)
            assert e.r_mean[i, k] == np.float32(stats.r_mean)
            assert e.r_dev[i, k] == np.float32(stats.r_dev)
        assert e.intensity_threshold[i] == np.float32(rec.intensity_threshold)
        assert e.gentrain_score[i] == np.float32(rec.cluster_score.total_score)
        assert e.cluster_sep[i] == np.float32(rec.cluster_score.cluster_separation)
        assert e.addresses[i] == rec.address


@pytest.mark.slow
def test_csv_manifest_consistent_with_bpm(gsa_bpm) -> None:  # type: ignore[no-untyped-def]  # IN-042
    need(GSA_CSV)
    c = read_csv_manifest(GSA_CSV)
    assert compare_manifests(gsa_bpm, c) == []
    # Assay type inferred from CSV sequences (fallback for BPM locus entries < v6, gtc2vcf rule) agrees with the
    # BPM's explicit value for 99.71% of Infinium I loci on GSA v4; the BPM is authoritative (known_unknowns K-003).
    inf1 = gsa_bpm.address_b != 0
    agree = (c.assay_types[inf1] == gsa_bpm.assay_types[inf1]).mean()
    assert agree >= 0.995


def test_gda_demo_sheet() -> None:
    need(GDA_SHEET)
    s = parse_samplesheet(GDA_SHEET)
    assert len(s.rows) == 36 and s.has_sample_name
    assert {r.values["Parent1"] for r in s.rows} >= {"NA12891"}