#!/usr/bin/env python3
"""Verify all complete SEQC2 WES BAM/BAI pairs with anonymous range GETs."""

from __future__ import annotations

import csv
import hashlib
import json
import statistics
import sys
from concurrent.futures import ThreadPoolExecutor, as_completed
from datetime import datetime, timezone
from pathlib import Path

from audit_public_wes_alignments import range_probe


DATE = "2026-08-01"
BASE = "https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/seqc/Somatic_Mutation_WG/data/WES/"


def read_md5(path: Path, prefix: str) -> dict[str, str]:
    out: dict[str, str] = {}
    for line in path.read_text(encoding="utf-8").splitlines():
        fields = line.split()
        if len(fields) < 2:
            continue
        filename = Path(fields[-1]).name
        if filename.startswith(prefix):
            out[filename] = fields[0]
    return out


def main() -> int:
    if len(sys.argv) != 2:
        raise SystemExit("usage: audit_seqc2_wes.py ANALYSIS_DIR")
    analysis = Path(sys.argv[1])
    bam_md5 = read_md5(analysis / f"seqc2_wes_bam_md5_{DATE}.txt", "WES_")
    bai_md5 = read_md5(analysis / f"seqc2_all_bai_md5_{DATE}.txt", "WES_")
    rows: list[dict[str, object]] = []
    for bam_name, md5 in sorted(bam_md5.items()):
        bai_name = bam_name[:-1] + "i"
        rows.append(
            {
                "repository": "NCBI_SEQC2",
                "collection": "HCC1395_HCC1395BL_WES",
                "sample_type": "tumor HCC1395" if "_T_" in bam_name else "matched normal HCC1395BL",
                "center_code": bam_name.split("_")[1],
                "replicate": bam_name.split("_")[3].split(".")[0],
                "reference_build": "GRCh38.d1.vd1",
                "capture": "Agilent SureSelect Human All Exon v6 + UTR",
                "alignment_url": BASE + bam_name,
                "index_url": BASE + bai_name,
                "alignment_md5": md5,
                "index_md5": bai_md5.get(bai_name, ""),
            }
        )

    tasks = {}
    # NCBI's public FTP gateway throttles bursts. Keep concurrency deliberately
    # low and allow a longer exponential retry window so 503s are not mistaken
    # for missing objects.
    with ThreadPoolExecutor(max_workers=2) as pool:
        for i, row in enumerate(rows):
            tasks[pool.submit(range_probe, str(row["alignment_url"]), None, 8)] = (i, "alignment")
            tasks[pool.submit(range_probe, str(row["index_url"]), None, 8)] = (i, "index")
        for future in as_completed(tasks):
            i, kind = tasks[future]
            rows[i][kind + "_probe"] = future.result()

    for row in rows:
        row["alignment_bytes"] = row["alignment_probe"].get("observed_total_bytes")  # type: ignore[index,union-attr]
        row["index_bytes"] = row["index_probe"].get("observed_total_bytes")  # type: ignore[index,union-attr]
        row["verified_downloadable"] = bool(
            row["alignment_probe"].get("ok") and row["index_probe"].get("ok")  # type: ignore[union-attr]
        )

    csv_path = analysis / f"seqc2_complete_wes_alignment_inventory_verified_{DATE}.csv"
    fields = [
        "repository", "collection", "sample_type", "center_code", "replicate",
        "reference_build", "capture", "alignment_url", "index_url", "alignment_bytes",
        "index_bytes", "alignment_md5", "index_md5", "alignment_http_status",
        "index_http_status", "verified_downloadable",
    ]
    with csv_path.open("w", encoding="utf-8", newline="") as fh:
        writer = csv.DictWriter(fh, fieldnames=fields)
        writer.writeheader()
        for row in rows:
            writer.writerow(
                {
                    **{k: row.get(k, "") for k in fields},
                    "alignment_http_status": row["alignment_probe"].get("status"),  # type: ignore[union-attr]
                    "index_http_status": row["index_probe"].get("status"),  # type: ignore[union-attr]
                }
            )

    good = [row for row in rows if row["verified_downloadable"]]
    sizes = [int(row["alignment_bytes"]) for row in rows if row["alignment_bytes"] is not None]
    payload = {
        "audit_date": DATE,
        "generated_at_utc": datetime.now(timezone.utc).isoformat(),
        "official_directory": BASE,
        "official_bam_md5_manifest": BASE + "WES.md5",
        "reference_paper": "https://pmc.ncbi.nlm.nih.gov/articles/PMC8532138/",
        "cnv_benchmark_paper": "https://pmc.ncbi.nlm.nih.gov/articles/PMC11188507/",
        "design": {
            "biological_materials": ["HCC1395 tumor", "HCC1395BL matched normal"],
            "alignment_objects": len(rows),
            "tumor_objects": sum(row["sample_type"] == "tumor HCC1395" for row in rows),
            "matched_normal_objects": sum(row["sample_type"] == "matched normal HCC1395BL" for row in rows),
            "capture": "Agilent SureSelect Human All Exon v6 + UTR",
            "reference_build": "GRCh38.d1.vd1",
            "caveat": "Technical multi-center replicates of two cell lines, not 24 independent patients; hyperdiploid somatic CNV benchmark, not germline exon-CNV truth.",
        },
        "verification": {
            "method": "anonymous HTTP GET Range bytes=0-0 on each BAM and BAI",
            "verified_pairs": len(good),
            "failed_pairs": len(rows) - len(good),
            "total_alignment_bytes": sum(sizes),
            "min_alignment_bytes": min(sizes),
            "median_alignment_bytes": int(statistics.median(sizes)),
            "max_alignment_bytes": max(sizes),
        },
        "matching_truth": {
            "description": "SEQC2 high-confidence HCC1395 somatic CNV call set integrated six callers across 21 WGS replicates and validated with Affymetrix CytoScan, Illumina BeadChip, and Bionano.",
            "fitness_for_LGASieve": "Strong real-data stress test for losses/gains, capture-site effects, replicate reproducibility, and tumor-normal normalization; not a direct clinical germline BRCA exon-LGA sensitivity cohort.",
        },
        "failures": [
            {
                "alignment_url": row["alignment_url"],
                "index_url": row["index_url"],
                "alignment_probe": row["alignment_probe"],
                "index_probe": row["index_probe"],
            }
            for row in rows
            if not row["verified_downloadable"]
        ],
        "csv": csv_path.name,
    }
    json_path = analysis / f"seqc2_complete_wes_repository_audit_{DATE}.json"
    json_path.write_text(json.dumps(payload, indent=2), encoding="utf-8")
    print(json.dumps(payload["verification"], indent=2))
    print(csv_path)
    print(json_path)
    return 0


if __name__ == "__main__":
    raise SystemExit(main())
