#!/usr/bin/env python3
"""Audit public complete-WES resources relevant to BRCA exon-CNV validation.

This intentionally distinguishes three concepts that are often conflated:

1. a complete WES alignment that is publicly downloadable;
2. raw complete WES reads that can be aligned, but are not a BAM/CRAM yet; and
3. a sample with independent BRCA1/2 exon-level deletion/duplication truth.

Every alignment and raw-read URL reported as downloadable is checked with both
HTTP HEAD and an actual 64-byte range GET.  No complete sequencing file is
downloaded by this audit.
"""

from __future__ import annotations

import base64
import csv
import datetime as dt
import hashlib
import html.parser
import json
import re
import time
import urllib.error
import urllib.parse
import urllib.request
from concurrent.futures import ThreadPoolExecutor, as_completed
from pathlib import Path
from typing import Any, Iterable


PROJECT_ROOT = Path(__file__).resolve().parents[1]
ANALYSIS_DIR = PROJECT_ROOT / "analysis"
ALIGNMENT_CSV = ANALYSIS_DIR / "public_full_wes_brca_specific_alignment_inventory_2026-08-01.csv"
RAW_CSV = ANALYSIS_DIR / "public_brca_candidate_raw_wes_inventory_2026-08-01.csv"
AUDIT_JSON = ANALYSIS_DIR / "public_full_wes_brca_specific_audit_2026-08-01.json"

USER_AGENT = "LGASieve-public-WES-audit/1.0"
TIMEOUT_SECONDS = 45
MAX_WORKERS = 16
RANGE_END = 63


def request_bytes(url: str, *, method: str = "GET", range_get: bool = False) -> dict[str, Any]:
    """HEAD or 64-byte GET with retries; return auditable response metadata."""

    headers = {"User-Agent": USER_AGENT}
    if range_get:
        headers["Range"] = f"bytes=0-{RANGE_END}"
    last_error = None
    for attempt in range(3):
        try:
            req = urllib.request.Request(url, headers=headers, method=method)
            with urllib.request.urlopen(req, timeout=TIMEOUT_SECONDS) as response:
                body = b"" if method == "HEAD" else response.read(RANGE_END + 1)
                return {
                    "status": response.status,
                    "content_length": _int_or_none(response.headers.get("Content-Length")),
                    "accept_ranges": response.headers.get("Accept-Ranges"),
                    "content_range": response.headers.get("Content-Range"),
                    "bytes_read": len(body),
                    "first_8_hex": body[:8].hex(),
                    "error": None,
                }
        except Exception as exc:  # recorded rather than hidden
            last_error = f"{type(exc).__name__}: {exc}"
            time.sleep(0.4 * (attempt + 1))
    return {
        "status": None,
        "content_length": None,
        "accept_ranges": None,
        "content_range": None,
        "bytes_read": 0,
        "first_8_hex": "",
        "error": last_error,
    }


def _int_or_none(value: Any) -> int | None:
    try:
        return int(value)
    except (TypeError, ValueError):
        return None


def fetch_json(url: str) -> Any:
    req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT})
    with urllib.request.urlopen(req, timeout=TIMEOUT_SECONDS) as response:
        return json.load(response)


def fetch_text(url: str) -> str:
    req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT})
    with urllib.request.urlopen(req, timeout=TIMEOUT_SECONDS) as response:
        return response.read().decode("utf-8", errors="replace")


def gcs_objects(prefix: str) -> list[dict[str, Any]]:
    base = "https://storage.googleapis.com/storage/v1/b/brain-genomics-public/o"
    token = None
    objects: list[dict[str, Any]] = []
    while True:
        params = {
            "prefix": prefix,
            "maxResults": "1000",
            "fields": "nextPageToken,items(name,size,md5Hash,crc32c,generation,updated)",
        }
        if token:
            params["pageToken"] = token
        page = fetch_json(base + "?" + urllib.parse.urlencode(params))
        objects.extend(page.get("items", []))
        token = page.get("nextPageToken")
        if not token:
            return objects


def gcs_url(name: str) -> str:
    return "https://storage.googleapis.com/brain-genomics-public/" + urllib.parse.quote(name, safe="/")


class LinkParser(html.parser.HTMLParser):
    def __init__(self) -> None:
        super().__init__()
        self.hrefs: list[str] = []

    def handle_starttag(self, tag: str, attrs: list[tuple[str, str | None]]) -> None:
        if tag.lower() != "a":
            return
        href = dict(attrs).get("href")
        if href:
            self.hrefs.append(href)


def list_directory_files(url: str) -> list[str]:
    parser = LinkParser()
    parser.feed(fetch_text(url))
    return [urllib.parse.urljoin(url, href) for href in parser.hrefs if not href.endswith("/")]


def sample_from_name(name: str) -> str:
    aliases = [
        (r"NA12878|HG001", "HG001/NA12878"),
        (r"NA24385|HG002", "HG002/NA24385"),
        (r"NA24149|HG003", "HG003/NA24149"),
        (r"NA24143|HG004", "HG004/NA24143"),
        (r"NA24631|HG005", "HG005/NA24631"),
        (r"HG006", "HG006/NA24694"),
        (r"HG007", "HG007/NA24695"),
        (r"NA12891", "NA12891"),
        (r"NA12892", "NA12892"),
        (r"NA18949", "NA18949"),
        (r"HG01528", "HG01528"),
    ]
    for pattern, sample in aliases:
        if re.search(pattern, name, re.IGNORECASE):
            return sample
    return "unknown"


def infer_gcs_metadata(path: str) -> dict[str, str]:
    parts = path.split("/")
    build = parts[2] if len(parts) > 2 else "unknown"
    platform = parts[4] if len(parts) > 4 else "unknown"
    capture = parts[5] if len(parts) > 5 else "unknown"
    depth = parts[6] if len(parts) > 6 else "unknown"
    return {
        "sample_id": sample_from_name(parts[-1]),
        "reference_build": build.upper(),
        "platform": platform,
        "capture": capture,
        "nominal_depth": depth,
    }


def brca_truth(sample: str) -> tuple[str, str, str, str]:
    if sample == "NA18949":
        return (
            "independent_positive",
            "BRCA1 heterozygous deletion; canonical exons E14-E15",
            "Phase 3 PASS SV + independent 30x ensemble PASS SV + DRAGEN Manta PASS SV",
            "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/integrated_sv_map/ALL.wgs.mergedSV.v8.20130502.svs.genotypes.vcf.gz",
        )
    if sample == "HG01528":
        return (
            "independent_positive",
            "BRCA1 heterozygous deletion; canonical exons E01-E06",
            "Independent 30x ensemble PASS SV + DRAGEN Manta PASS SV + DRAGEN depth CNV PASS CN=1",
            "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20210124.SV_Illumina_Integration/1KGP_3202.gatksv_svtools_novelins.freeze_V3.wAF.vcf.gz",
        )
    return (
        "not_gold_labeled_for_BRCA_exon_CNV",
        "",
        "No independent exon-level BRCA deletion/duplication truth located. For HG002, the GIAB v0.6 and CMRG v1.00 SV benchmark BEDs do not cover BRCA1 or BRCA2, so absence of a benchmark call is not a validated negative.",
        "https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/CMRG_v1.00/",
    )


def google_brain_alignment_pairs() -> list[dict[str, Any]]:
    rows: list[dict[str, Any]] = []
    for build in ("grch37", "grch38"):
        objects = gcs_objects(f"research/sequencing/{build}/bam/")
        selected = {
            obj["name"]: obj
            for obj in objects
            if re.search(r"/(?:wes|exome)", obj["name"], re.IGNORECASE)
            and (obj["name"].endswith(".bam") or obj["name"].endswith(".bam.bai"))
        }
        for name, bam in sorted(selected.items()):
            if not name.endswith(".bam"):
                continue
            bai = selected.get(name + ".bai")
            meta = infer_gcs_metadata(name)
            truth_status, event, evidence, evidence_url = brca_truth(meta["sample_id"])
            rows.append(
                {
                    "repository": "Google Brain Genomics Public / GIAB",
                    "dataset": "Baid et al. multi-kit gold-standard sequencing dataset",
                    **meta,
                    "bam_url": gcs_url(name),
                    "bam_bytes": int(bam["size"]),
                    "bam_md5_base64": bam.get("md5Hash", ""),
                    "bai_url": gcs_url(name + ".bai") if bai else "",
                    "bai_bytes": int(bai["size"]) if bai else None,
                    "bai_md5_base64": bai.get("md5Hash", "") if bai else "",
                    "brca_truth_status": truth_status,
                    "brca_truth_event": event,
                    "brca_truth_evidence": evidence,
                    "truth_evidence_url": evidence_url,
                    "validation_use": "multi-kit/depth/platform robustness only; not a BRCA exon-CNV truth control",
                }
            )
    return rows


def ncbi_giab_alignment_pairs() -> list[dict[str, Any]]:
    roots = [
        ("Oslo Agilent SureSelect v5", "GRCh37/hs37d5", "Illumina HiSeq", "100-200x", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG002_NA24385_son/OsloUniversityHospital_Exome/"),
        ("Oslo Agilent SureSelect v5", "GRCh37/hs37d5", "Illumina HiSeq", "100-200x", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG003_NA24149_father/OsloUniversityHospital_Exome/"),
        ("Oslo Agilent SureSelect v5", "GRCh37/hs37d5", "Illumina HiSeq", "100-200x", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG004_NA24143_mother/OsloUniversityHospital_Exome/"),
        ("Oslo Agilent SureSelect v5", "GRCh37/hs37d5", "Illumina HiSeq", "100-200x", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/ChineseTrio/HG005_NA24631_son/OsloUniversityHospital_Exome/"),
        ("Ion AmpliSeq Exome", "hg19", "Ion Proton", "mean 256x per run", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG002_NA24385_son/ion_exome/"),
        ("Ion AmpliSeq Exome", "hg19", "Ion Proton", "mean 256x per run", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG003_NA24149_father/ion_exome/"),
        ("Ion AmpliSeq Exome", "hg19", "Ion Proton", "mean 256x per run", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG004_NA24143_mother/ion_exome/"),
        ("Ion AmpliSeq Exome", "hg19", "Ion Proton", "mean 256x per run", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/ChineseTrio/HG005_NA24631_son/ion_exome/"),
        ("Garvan Nextera Rapid Capture Exome", "hg19", "Illumina HiSeq 2500", "per-lane", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/NA12878/Garvan_NA12878_HG001_HiSeq_Exome/"),
        ("Nebraska TruSeq Exome", "hg19", "Illumina", "unspecified", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/NA12878/Nebraska_NA12878_HG001_TruSeq_Exome/"),
        ("Ion exome", "b37 or hg19 (filename-specific)", "Ion Torrent", "unspecified", "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/NA12878/ion_exome/"),
    ]
    rows: list[dict[str, Any]] = []
    for capture, build, platform, depth, root in roots:
        files = set(list_directory_files(root))
        for bam_url in sorted(url for url in files if url.endswith(".bam")):
            candidates = [bam_url + ".bai", bam_url[:-4] + ".bai"]
            bai_url = next((url for url in candidates if url in files), "")
            sample = sample_from_name(bam_url)
            truth_status, event, evidence, evidence_url = brca_truth(sample)
            if capture == "Ion exome":
                this_build = "b37" if ".b37." in bam_url else "hg19"
            else:
                this_build = build
            rows.append(
                {
                    "repository": "NIST/NCBI Genome in a Bottle",
                    "dataset": "GIAB legacy complete exome alignments",
                    "sample_id": sample,
                    "reference_build": this_build,
                    "platform": platform,
                    "capture": capture,
                    "nominal_depth": depth,
                    "bam_url": bam_url,
                    "bam_bytes": None,
                    "bam_md5_base64": "",
                    "bai_url": bai_url,
                    "bai_bytes": None,
                    "bai_md5_base64": "",
                    "brca_truth_status": truth_status,
                    "brca_truth_event": event,
                    "brca_truth_evidence": evidence,
                    "truth_evidence_url": evidence_url,
                    "validation_use": "technical WES control only; not a BRCA exon-CNV truth control",
                }
            )
    return rows


def positive_1000g_alignment_pairs() -> list[dict[str, Any]]:
    base = "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data"
    specs = [
        ("NA18949", "JPT", "NA18949.mapped.ILLUMINA.bwa.JPT.exome.20121211.bam", 11042427115, 7816264, "107.282x"),
        ("HG01528", "IBS", "HG01528.mapped.ILLUMINA.bwa.IBS.exome.20121211.bam", 10319084612, 6814936, "138.939x"),
    ]
    rows = []
    for sample, population, filename, bam_bytes, bai_bytes, depth in specs:
        bam_url = f"{base}/{sample}/exome_alignment/{filename}"
        truth_status, event, evidence, evidence_url = brca_truth(sample)
        rows.append(
            {
                "repository": "1000 Genomes Project Phase 3",
                "dataset": "official Phase 3 complete exome alignment",
                "sample_id": sample,
                "reference_build": "GRCh37",
                "platform": "Illumina",
                "capture": f"1000G Phase 3 exome ({population})",
                "nominal_depth": depth,
                "bam_url": bam_url,
                "bam_bytes": bam_bytes,
                "bam_md5_base64": "",
                "bai_url": bam_url + ".bai",
                "bai_bytes": bai_bytes,
                "bai_md5_base64": "",
                "brca_truth_status": truth_status,
                "brca_truth_event": event,
                "brca_truth_evidence": evidence,
                "truth_evidence_url": evidence_url,
                "validation_use": "independent BRCA1 exon-deletion positive control",
            }
        )
    return rows


def audit_urls(urls: Iterable[str]) -> dict[str, dict[str, Any]]:
    unique = sorted({url for url in urls if url})

    def one(url: str) -> tuple[str, dict[str, Any]]:
        head = request_bytes(url, method="HEAD")
        ranged = request_bytes(url, method="GET", range_get=True)
        return url, {"head": head, "range_get": ranged}

    # NCBI's FTP-over-HTTPS front end rate-limits parallel range requests and
    # emits transient 503s.  Keep that host deliberately serial while retaining
    # bounded concurrency for object stores designed for it.
    serial_ncbi = [url for url in unique if urllib.parse.urlparse(url).hostname == "ftp.ncbi.nlm.nih.gov"]
    parallel = [url for url in unique if url not in set(serial_ncbi)]
    results: dict[str, dict[str, Any]] = {}

    def run_group(group: list[str], workers: int) -> None:
        with ThreadPoolExecutor(max_workers=workers) as pool:
            futures = {pool.submit(one, url): url for url in group}
            for future in as_completed(futures):
                url = futures[future]
                try:
                    key, value = future.result()
                except Exception as exc:
                    key, value = url, {"head": {"error": str(exc)}, "range_get": {"error": str(exc)}}
                results[key] = value

    run_group(parallel, MAX_WORKERS)
    run_group(serial_ncbi, 2)
    return results


def attach_audit(rows: list[dict[str, Any]], audits: dict[str, dict[str, Any]]) -> None:
    for row in rows:
        bam_audit = audits.get(row["bam_url"], {})
        bai_audit = audits.get(row["bai_url"], {}) if row.get("bai_url") else {}
        bh = bam_audit.get("head", {})
        br = bam_audit.get("range_get", {})
        ih = bai_audit.get("head", {})
        ir = bai_audit.get("range_get", {})
        if row.get("bam_bytes") is None:
            row["bam_bytes"] = bh.get("content_length")
        if row.get("bai_bytes") is None:
            row["bai_bytes"] = ih.get("content_length")
        row.update(
            {
                "bam_head_status": bh.get("status"),
                "bam_head_length": bh.get("content_length"),
                "bam_accept_ranges": bh.get("accept_ranges"),
                "bam_range_status": br.get("status"),
                "bam_content_range": br.get("content_range"),
                "bam_range_bytes_read": br.get("bytes_read"),
                "bam_first_8_hex": br.get("first_8_hex"),
                "bam_error": bh.get("error") or br.get("error"),
                "bai_head_status": ih.get("status"),
                "bai_head_length": ih.get("content_length"),
                "bai_accept_ranges": ih.get("accept_ranges"),
                "bai_range_status": ir.get("status"),
                "bai_content_range": ir.get("content_range"),
                "bai_range_bytes_read": ir.get("bytes_read"),
                "bai_first_8_hex": ir.get("first_8_hex"),
                "bai_error": ih.get("error") or ir.get("error"),
            }
        )
        expected_bam = row.get("bam_bytes")
        expected_bai = row.get("bai_bytes")
        row["pair_downloadable"] = bool(
            bh.get("status") == 200
            and br.get("status") == 206
            and br.get("bytes_read") == RANGE_END + 1
            and row.get("bai_url")
            and ih.get("status") == 200
            and ir.get("status") == 206
            and ir.get("bytes_read") == RANGE_END + 1
            and (expected_bam is None or bh.get("content_length") == expected_bam)
            and (expected_bai is None or ih.get("content_length") == expected_bai)
        )


BRACNAC_CANDIDATES = {
    "SRR5604273": "del_BRCA1_ex13-ex13; authors suggest possible germline CNV",
    "SRR5604275": "del_BRCA1_ex2-ex23; likely somatic CNA",
    "SRR5604279": "del_BRCA1_ex17-ex21 / possibly whole-gene loss; likely somatic CNA",
    "SRR5604281": "dupl_BRCA1_ex13-ex23; likely somatic CNA",
    "SRR5604292": "dupl_BRCA1_ex2-ex23; likely somatic CNA",
    "SRR5604295": "del_BRCA1_ex13-ex19; authors suggest possible germline CNV",
    "SRR5604298": "del_BRCA1_ex11-ex12; likely somatic CNA",
    "SRR5604299": "del_BRCA1_ex2-ex23; likely somatic CNA",
    "SRR5604308": "dupl_BRCA1_ex17-ex22; likely somatic CNA",
    "SRR5604312": "del_BRCA2_ex15-ex18; likely somatic CNA",
    "SRR5604313": "dupl_BRCA1_ex15-ex16; likely somatic CNA",
    "SRR5604314": "dupl_BRCA1_ex13-ex17 and ex19-ex23; likely somatic CNA",
    "SRR5604315": "del_BRCA1_ex12-ex13; likely somatic CNA",
}


def ena_raw_wes_rows() -> list[dict[str, Any]]:
    fields = ",".join(
        [
            "run_accession",
            "experiment_accession",
            "sample_accession",
            "scientific_name",
            "library_strategy",
            "library_source",
            "library_layout",
            "instrument_platform",
            "fastq_ftp",
            "fastq_bytes",
            "fastq_md5",
            "submitted_ftp",
            "submitted_bytes",
            "submitted_md5",
        ]
    )
    url = (
        "https://www.ebi.ac.uk/ena/portal/api/filereport?"
        + urllib.parse.urlencode(
            {
                "accession": "PRJNA388048",
                "result": "read_run",
                "fields": fields,
                "format": "json",
            }
        )
    )
    runs = fetch_json(url)
    rows: list[dict[str, Any]] = []
    for run in runs:
        urls = run.get("fastq_ftp", "").split(";") if run.get("fastq_ftp") else []
        sizes = run.get("fastq_bytes", "").split(";") if run.get("fastq_bytes") else []
        md5s = run.get("fastq_md5", "").split(";") if run.get("fastq_md5") else []
        for ordinal, fq_url in enumerate(urls, start=1):
            run_id = run["run_accession"]
            rows.append(
                {
                    "repository": "ENA/NCBI SRA",
                    "bioproject": "PRJNA388048",
                    "run_accession": run_id,
                    "experiment_accession": run.get("experiment_accession", ""),
                    "sample_accession": run.get("sample_accession", ""),
                    "library_strategy": run.get("library_strategy", ""),
                    "library_layout": run.get("library_layout", ""),
                    "file_ordinal": ordinal,
                    "file_url": "https://" + fq_url,
                    "file_bytes": int(sizes[ordinal - 1]) if ordinal <= len(sizes) and sizes[ordinal - 1] else None,
                    "file_md5": md5s[ordinal - 1] if ordinal <= len(md5s) else "",
                    "paper_candidate": run_id in BRACNAC_CANDIDATES,
                    "paper_candidate_event": BRACNAC_CANDIDATES.get(run_id, ""),
                    "truth_status": "computational_candidate_not_independently_validated" if run_id in BRACNAC_CANDIDATES else "unlabeled",
                    "alignment_status": "raw paired FASTQ only; no public submitted BAM/CRAM in ENA file report",
                }
            )
    return rows


def attach_raw_audit(rows: list[dict[str, Any]], audits: dict[str, dict[str, Any]]) -> None:
    for row in rows:
        audit = audits.get(row["file_url"], {})
        head = audit.get("head", {})
        ranged = audit.get("range_get", {})
        row.update(
            {
                "head_status": head.get("status"),
                "head_length": head.get("content_length"),
                "accept_ranges": head.get("accept_ranges"),
                "range_status": ranged.get("status"),
                "content_range": ranged.get("content_range"),
                "range_bytes_read": ranged.get("bytes_read"),
                "error": head.get("error") or ranged.get("error"),
            }
        )
        row["downloadable"] = bool(
            head.get("status") == 200
            and ranged.get("status") == 206
            and ranged.get("bytes_read") == RANGE_END + 1
            and head.get("content_length") == row.get("file_bytes")
        )


def write_csv(path: Path, rows: list[dict[str, Any]]) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    fields: list[str] = []
    seen: set[str] = set()
    for row in rows:
        for key in row:
            if key not in seen:
                fields.append(key)
                seen.add(key)
    with path.open("w", encoding="utf-8", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=fields)
        writer.writeheader()
        writer.writerows(rows)


def sha256(path: Path) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as handle:
        for chunk in iter(lambda: handle.read(1024 * 1024), b""):
            digest.update(chunk)
    return digest.hexdigest()


def main() -> None:
    alignment_rows = (
        positive_1000g_alignment_pairs()
        + ncbi_giab_alignment_pairs()
        + google_brain_alignment_pairs()
    )
    alignment_audits = audit_urls(
        url
        for row in alignment_rows
        for url in (row.get("bam_url", ""), row.get("bai_url", ""))
        if url
    )
    attach_audit(alignment_rows, alignment_audits)
    write_csv(ALIGNMENT_CSV, alignment_rows)

    raw_rows = ena_raw_wes_rows()
    raw_audits = audit_urls(row["file_url"] for row in raw_rows)
    attach_raw_audit(raw_rows, raw_audits)
    write_csv(RAW_CSV, raw_rows)

    positive_rows = [row for row in alignment_rows if row["brca_truth_status"] == "independent_positive"]
    google_rows = [row for row in alignment_rows if row["repository"].startswith("Google Brain")]
    nist_rows = [row for row in alignment_rows if row["repository"].startswith("NIST/NCBI")]
    summary = {
        "schema_version": "1.0",
        "generated_utc": dt.datetime.now(dt.timezone.utc).isoformat(),
        "scope": "Public complete human WES BAM/BAI pairs with special emphasis on independently supported BRCA1/2 exon-level deletion/duplication truth; plus a separately labelled raw-WES candidate cohort.",
        "methods": {
            "downloadability_definition": "HEAD HTTP 200 with exact expected Content-Length plus actual GET Range bytes=0-63 returning HTTP 206 and exactly 64 bytes, for both BAM and BAI.",
            "full_file_downloaded": False,
            "range_bytes_per_file": RANGE_END + 1,
            "alignment_pairs_audited": len(alignment_rows),
            "alignment_urls_audited": len(alignment_audits),
            "raw_fastq_files_audited": len(raw_rows),
        },
        "results": {
            "complete_wes_bam_pairs_found_in_this_brca_specific_search": len(alignment_rows),
            "complete_wes_bam_pairs_downloadable_now": sum(bool(row["pair_downloadable"]) for row in alignment_rows),
            "complete_wes_bam_pairs_failed_audit": sum(not bool(row["pair_downloadable"]) for row in alignment_rows),
            "unique_complete_wes_bam_content_md5_where_available": len({row["bam_md5_base64"] for row in alignment_rows if row.get("bam_md5_base64")}),
            "exact_global_unique_bam_content_count": None,
            "exact_global_unique_bam_content_count_reason": "Only the 266 Google Brain objects expose a uniform content MD5 in the enumerated metadata. Their 242 unique MD5s are exact within that repository; the 28 NIST and 2 1000G objects are kept separate rather than guessed to be globally unique.",
            "distinct_people_in_alignment_inventory": len({row["sample_id"] for row in alignment_rows}),
            "independently_supported_BRCA_exon_CNV_positive_people": len({row["sample_id"] for row in positive_rows}),
            "independently_supported_BRCA_exon_CNV_positive_alignments": len(positive_rows),
            "confirmed_deletion_positive_people": 2,
            "confirmed_duplication_positive_people": 0,
            "google_brain_giab_bams": len(google_rows),
            "google_brain_giab_bams_grch37": sum(row["reference_build"] == "GRCH37" for row in google_rows),
            "google_brain_giab_bams_grch38": sum(row["reference_build"] == "GRCH38" for row in google_rows),
            "google_brain_unique_bam_content_md5": len({row["bam_md5_base64"] for row in google_rows if row.get("bam_md5_base64")}),
            "google_brain_total_bam_bytes": sum(row["bam_bytes"] or 0 for row in google_rows),
            "nist_giab_legacy_bams": len(nist_rows),
            "raw_PRJNA388048_runs": len({row["run_accession"] for row in raw_rows}),
            "raw_PRJNA388048_fastq_files": len(raw_rows),
            "raw_PRJNA388048_fastq_files_downloadable_now": sum(bool(row["downloadable"]) for row in raw_rows),
            "raw_PRJNA388048_total_fastq_bytes": sum(row["file_bytes"] or 0 for row in raw_rows),
            "raw_PRJNA388048_BRACNAC_candidate_runs": len({row["run_accession"] for row in raw_rows if row["paper_candidate"]}),
        },
        "positive_controls": [
            {
                "sample_id": row["sample_id"],
                "event": row["brca_truth_event"],
                "bam_url": row["bam_url"],
                "bam_bytes": row["bam_bytes"],
                "bai_url": row["bai_url"],
                "bai_bytes": row["bai_bytes"],
                "pair_downloadable": row["pair_downloadable"],
                "truth_evidence": row["brca_truth_evidence"],
                "truth_evidence_url": row["truth_evidence_url"],
            }
            for row in positive_rows
        ],
        "research_candidate_not_truth": {
            "dataset": "PRJNA388048",
            "description": "60 public tumor WES runs; BRACNAC reported 13 BRCA CNA candidate runs, but the WES validation column was explicitly 'No' and the calls were not independently/orthogonally validated. ENA exposes paired FASTQs, not ready-to-query indexed BAM/CRAM files.",
            "candidate_runs": BRACNAC_CANDIDATES,
            "paper": "https://pmc.ncbi.nlm.nih.gov/articles/PMC10706169/",
            "bioproject": "https://www.ncbi.nlm.nih.gov/bioproject/388048",
        },
        "searched_but_not_qualifying": [
            {
                "resource": "TCGA-BRCA and TCGA-OV",
                "reason": "Individual WES BAMs are controlled/protected and require dbGaP authorization; public CNV matrices/VCFs are not complete WES read files.",
                "evidence": "https://gdc.cancer.gov/access-data/data-access-processes-and-tools",
            },
            {
                "resource": "EGA BRCA-mutant tumor WES cohorts (including EGAD00001008276)",
                "reason": "Controlled access through a Data Access Committee; cannot be downloaded anonymously now.",
                "evidence": "https://ega-archive.org/datasets/EGAD00001008276",
            },
            {
                "resource": "HSCGP WES (phs001343 / EGAS00001002400)",
                "reason": "Portal explicitly states WES primary data are controlled access.",
                "evidence": "https://hscgp.broadinstitute.org/",
            },
            {
                "resource": "Northern Finnish hereditary breast cancer WES (98 cases)",
                "reason": "Paper explicitly says individual exome data cannot be made publicly available due to local regulation.",
                "evidence": "https://doi.org/10.1371/journal.pgen.1010889",
            },
            {
                "resource": "Horizon BRCA/reference-standard exome downloads",
                "reason": "The public downloads are VCFs, not complete WES BAM/CRAM or FASTQ reads.",
                "evidence": "https://horizondiscovery.com/en/resources/exome-data",
            },
            {
                "resource": "Seraseq BRCA1/2 LGR reference materials",
                "reason": "Commercial material; the LGR constructs are plasmid spikes and the manufacturer states they are incompatible with coverage-depth CNV analysis because they do not create genome-wide copy loss/gain.",
                "evidence": "https://www.seracare.com/globalassets/seracare-resources/pi-seraseq-ffpe-brca1-2-lgr-ref-mtrl.pdf",
            },
            {
                "resource": "gnomAD / ExAC / UK Biobank / GeneDx aggregate resources",
                "reason": "Aggregate variants or controlled individual data; no anonymously downloadable individual complete WES BAM/CRAM paired with BRCA exon-CNV truth.",
                "evidence": "https://gnomad.broadinstitute.org/help/terms-of-use",
            },
            {
                "resource": "ClinVar / BRCA Exchange / dbVar",
                "reason": "Variant assertions/call records do not provide matching complete individual WES alignments and indexes.",
                "evidence": "https://www.ncbi.nlm.nih.gov/clinvar/",
            },
        ],
        "limitations": [
            "The inventory is an extensive reproducible search, not a proof that no qualifying file exists anywhere on the internet.",
            "A downloadable GIAB exome is not automatically a BRCA CNV negative. The HG002 SV benchmark confidence intervals examined here exclude BRCA1 and BRCA2.",
            "The 1000 Genomes/DRAGEN truth is research-grade computational corroboration, not a clinical MLPA/ddPCR truth set.",
            "Only two independently supported deletion-positive people and no duplication-positive person were found with ready-to-download complete WES BAM/BAI files.",
            "Downsampled or alternate-reference BAMs from the same biological person are technical replicates, not independent clinical samples.",
        ],
        "artifacts": {
            "alignment_inventory_csv": ALIGNMENT_CSV.name,
            "raw_candidate_inventory_csv": RAW_CSV.name,
        },
    }
    AUDIT_JSON.write_text(json.dumps(summary, indent=2, ensure_ascii=False) + "\n", encoding="utf-8")
    summary["artifacts"]["alignment_inventory_sha256"] = sha256(ALIGNMENT_CSV)
    summary["artifacts"]["raw_candidate_inventory_sha256"] = sha256(RAW_CSV)
    AUDIT_JSON.write_text(json.dumps(summary, indent=2, ensure_ascii=False) + "\n", encoding="utf-8")
    print(json.dumps(summary["results"], indent=2))
    print(f"alignment_csv={ALIGNMENT_CSV}")
    print(f"raw_csv={RAW_CSV}")
    print(f"audit_json={AUDIT_JSON}")


if __name__ == "__main__":
    main()
