#!/usr/bin/env python3
"""collect_summary.py: read one sample's step outputs into one JSON summary.

Both reports (scripts/24-html-report.sh and scripts/generate-report.sh) are
rendered from this summary by bin/render_report.py, so they cannot disagree
about a count. Standard library only (it runs in the plain python image).

  collect_summary.py --sample S --sample-dir DIR [--out S/summary.json]

DIR is a bash output folder (GENOME_DIR/S) or a Nextflow one (outdir/S): both
publish names are read (vcf/ or roh/ for ROH, mosdepth/ or coverage/ for
depth, hla_t1k/ or hla/ for HLA, vcf/ or pharmcat/ for PharmCAT).

Each section has a state:
  ok          read from its file
  missing     no file: the step has not run
  failed      no file, and its step failed in the latest run-all.sh run
              (logs/run_status.tsv)
  stale       the file is older than the latest run-all.sh run and its step was
              not ok in that run (logs/run_status.tsv); the values are shown
              with the file's date
  unreadable  the file exists but could not be read (the reason is in note)

The run manifest (S/run_manifest.tsv, bin/write_manifest.sh) is copied in for
the report footer, with the ClinVar release date and the HLA database release.
Two optional keys, for the report's run facts and badges: recheck, the last
database re-check (S/recheck/state.json, bin/owngenome_recheck.py), and new_since,
the date of the previous run kept in S/history/ and the sections whose results
changed since it (bin/run_delta.py's comparison).

  collect_summary.py sample-qc --sample S --somalier-samples F --selfsm F [...]

writes the sample identity and contamination table of step 33 (and of the
Nextflow SAMPLE_QC process): the sex somalier infers from the reads, the
declared sex, VerifyBamID2's FREEMIX against a warning threshold, and the
other samples somalier finds to be the same person. Both callers stop on the
sex_check value this writes, so the rule lives here once.

  collect_summary.py prs-format --scores DIR --labels assets/pgs_scores.tsv --out DIR --alleles F
                                [--keep-x] [--keep-ambiguous --reference FASTA]
  collect_summary.py prs-table --sample S --results DIR --sampleset NAME --scores DIR --input-kind gvcf --out F
                               [--sex male|female]

are the two ends of step 25 (and of the Nextflow PRS_PREPARE and PRS_SUMMARY
processes) around pgsc_calc: prs-format writes each GRCh38-harmonised PGS
Catalog file as the custom GRCh38 scoring file pgsc_calc reads, labelled with
the catalog's trait, and lists every effect and other allele so the score
positions can be genotyped from the gVCF; prs-table reads pgsc_calc's match
summary and scores into <sample>_prs_summary.tsv, with the ancestry-adjusted
percentile and the reference group when pgsc_calc ran with an ancestry panel,
and then writes the principal components and population of step 26.

  collect_summary.py mito-deletions --eklipse eKLIPse_deletions.csv --out S_chrM_deletions.tsv
  collect_summary.py mito-copies --idxstats S_idxstats.txt --out S_mtdna_copies.tsv

are step 20's large mtDNA deletions (eKLIPse's CSV as a table) and its mtDNA
copy number (chrM against autosomal reads per base, from samtools idxstats),
for scripts/20-mtoolbox.sh and the Nextflow MITO_DELETIONS process.
"""
import argparse
import base64
import csv
import glob
import gzip
import json
import math
import os
import re
import struct
import sys
import zlib
from datetime import datetime, timezone

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import pgx_outside_calls  # noqa: E402  step 36's rule for a CYP2D6 call and an HLA allele
import ancestry_fine_report  # noqa: E402  step 26b's group names and painting colours
import paraphase_report  # noqa: E402  step 35b's tables and its SMN1 rule beside Parascopy
import clip_check  # noqa: E402  CPSR's variant ids and class names, as step 17's read check reads them
import annotsv_postfilter  # noqa: E402  step 05's class 4-5 checks and the step 22 support join
import chip_report  # noqa: E402  step 49's detection limit at the median depth
import run_delta  # noqa: E402  the comparison with the previous run, for the "New since" badges
import traits_report  # noqa: E402  step 42's genotype-rules columns and the rules the traits table holds
import genotype_rules_indels  # noqa: E402  the APOL1 G2 deletion beside step 42's APOL1 rule

# 2: cyp2d6.consensus, Cyrius's Filter applied, HLA quality per allele,
# Stranger's status on the repeat loci, the CPIC unclassified count, and the
# section state 'failed'.
SCHEMA_VERSION = 2

# ClinVar review status -> stars. bin/clinvar_hits.awk holds the same table
# (tests/test_collect_summary.py checks that the two agree).
STARS = {
    "practice_guideline": 4,
    "reviewed_by_expert_panel": 3,
    "criteria_provided,_multiple_submitters,_no_conflicts": 2,
    "criteria_provided,_single_submitter": 1,
    "criteria_provided,_conflicting_classifications": 1,
    "criteria_provided,_conflicting_interpretations": 1,
}

# The call floor of step 6 (scripts/06-clinvar-screen.sh writes a floor other
# than this default into the hits VCF as ##PGP_CALL_FLOOR=GQ=..,DP=..): a P/LP
# hit below either is a low-confidence call, listed but counted apart. A hits
# file without the line used this default, or was written before the floor.
CALL_FLOOR = {"GQ": 30, "DP": 10}
# gnomAD allele frequency above which a P/LP hit is marked common.
COMMON_AF = 0.01
# The genotype check of the clinical filter (bin/csq_rank.awk, step 23 and the
# CLINICAL_FILTER module): step 6's GQ and DP floor, and the alt read fraction
# (AD) a heterozygote must lie within and a homozygous alt call must reach.
# bin/csq_rank.awk writes the floors it used into the clinical VCF as
# ##PGP_CALL_FLOOR=GQ=..,DP=..,HET_AB=lo-hi,HOM_AB=..; these are its defaults.
AB_FLOOR = {"HET_AB": (0.2, 0.8), "HOM_AB": 0.8}
IMPACT_RANK = {"HIGH": 4, "MODERATE": 3, "LOW": 2, "MODIFIER": 1}

# Steps run-all.sh runs and records in logs/run_status.tsv: the rows of
# assets/steps.tsv with recorded=yes (tests/test_run_state.py keeps the two
# equal; the set stays written here because this file also runs where only
# bin/ is mounted). A section whose step is not here (variant calling, step
# 03, which run-all reuses on purpose) is never marked stale.
RUN_ALL_STEPS = {"03a", "03b", "03c", "03d", "04", "04b", "05", "06", "07", "08", "09", "09b", "10", "11",
                 "12", "13", "14", "15", "16", "16b", "17", "18", "19", "20", "21", "22", "23", "25", "26",
                 "26b", "27", "28", "29", "30", "31", "32", "33", "35", "35b", "37", "38", "39", "40", "41",
                 "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53"}

# What this pipeline does not assess, whatever ran.
NOT_ASSESSED = [
    "Mosaic and low-fraction small variants (only chrM heteroplasmy is reported, and with step 49 the CHIP genes)",
    "Methylation and imprinting",
    "Repeat expansions outside the ExpansionHunter catalog",
    "Phase: two variants in one gene are compound-het candidates, not confirmed",
    "Variants of uncertain significance (not interpreted)",
]

# Added to NOT_ASSESSED unless step 35's depth count (bin/hba_copy_number.py,
# short reads) or step 35b (Paraphase, PacBio HiFi, both genes) made a call.
HBA_NOT_ASSESSED = ("HBA1/HBA2 copy number (alpha-thalassaemia deletions): neither step 35's read-depth count "
                    "(short reads) nor step 35b (Paraphase, PacBio HiFi) made a call")

# Added to NOT_ASSESSED when step 50 (opt-in, MEGAnE) did not run; a failed
# step 50 is listed under not run, and unreadable tables under their section.
MEI_NOT_ASSESSED = ("Mobile element insertions (Alu, L1 and SVA, the FKTN, TAF1 and BBS1 SVA insertions among them): "
                    "not looked for (step 50, MEGAnE, is opt-in)")

# Added to NOT_ASSESSED unless step 51 (MoChA) assessed the sample.
MCA_NOT_ASSESSED = ("Autosomal mosaic chromosomal alterations (large gains, losses and copy-neutral LOH in a share "
                    "of the blood cells): step 51 (MoChA) did not assess this sample")

# Added to NOT_ASSESSED unless step 52 (VNtyper) gave a positive, possible or
# negative reading: no other step can see a frameshift inside the MUC1 VNTR.
VNTYPER_NOT_ASSESSED = ("MUC1 VNTR frameshifts (ADTKD-MUC1): step 52 (VNtyper, opt-in) gave no reading; "
                        "no standard caller sees them")

# Added to NOT_ASSESSED unless step 25 ran with an ancestry reference panel.
PRS_NOT_ADJUSTED = "Polygenic score percentiles: no ancestry reference panel was used, so the scores are raw sums"

EH_LOCI = ["HTT", "FMR1", "C9ORF72", "ATXN1", "DMPK"]

# ACMG SF v3.3 (Lee et al., Genet Med 2025;27(8):101454): the 84 genes in which
# the ACMG recommends reporting pathogenic and likely pathogenic variants found
# by chance, and the reporting rule of each gene (Table 1), are the rows of
# assets/acmg_sf.tsv. A new list version means new rows there and a new
# ACMG_SF_VERSION. The sample's carrier list reads assets/carrier_genes.tsv.
# Both are read from ASSETS_DIR: the repository's assets/ folder beside bin/,
# or the folder render_report.py --assets names (the containers of step 24 and
# of HTML_REPORT mount bin/ only).
ACMG_SF_VERSION = "ACMG SF v3.3"
ASSETS_DIR = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "assets")
SF_REPORT_IF = ("any", "biallelic", "hom_variant", "hemizygous", "het_female_too")
SEVERITY_KINDS = ("hypomorph", "pseudodeficiency")


class Unreadable(Exception):
    pass


def iso(ts):
    return datetime.fromtimestamp(ts, tz=timezone.utc).strftime("%Y-%m-%d %H:%M UTC")


def open_text(path):
    if path.endswith((".gz", ".bgz")):
        return gzip.open(path, "rt", errors="replace")
    return open(path, errors="replace")


def first_existing(base, rels):
    for rel in rels:
        for p in sorted(glob.glob(os.path.join(base, rel), recursive=True)):
            if os.path.isfile(p):
                return p
    return None


def newest(base, rels):
    found = [p for rel in rels for p in glob.glob(os.path.join(base, rel), recursive=True) if os.path.isfile(p)]
    return max(found, key=os.path.getmtime) if found else None


def vcf_records(path):
    """Yield the split data lines of a VCF (plain or gzip)."""
    try:
        with open_text(path) as f:
            for line in f:
                if line.startswith("#"):
                    continue
                yield line.rstrip("\n").split("\t")
    except (OSError, EOFError, gzip.BadGzipFile, zlib.error) as e:
        raise Unreadable(f"{os.path.basename(path)}: {e}") from e


def info_map(field):
    out = {}
    for kv in field.split(";"):
        k, _, v = kv.partition("=")
        out[k] = v
    return out


def read_tsv(path):
    with open_text(path) as f:
        return list(csv.DictReader(f, delimiter="\t"))


# --- sections -------------------------------------------------------------------

def sec_variants(d, s):
    p = first_existing(d, [f"vcf/{s}.vcf.gz"])
    if not p:
        return None, {}
    total = passed = snps = indels = 0
    for r in vcf_records(p):
        if len(r) < 7:
            continue
        total += 1
        if r[6] == "PASS":
            passed += 1
        alts = [a for a in r[4].split(",") if a not in (".", "*", "<NON_REF>")]
        if any(len(a) == len(r[3]) == 1 for a in alts):
            snps += 1
        if any(len(a) != len(r[3]) and not a.startswith("<") for a in alts):
            indels += 1
    return p, {"total": total, "pass": passed, "snps": snps, "indels": indels}


def clinvar_variants(v, sep):
    """'chr1:10:C:T|chr1:15:G:A' -> 'chr1:10 C>T<sep>chr1:15 G>A'."""
    out = []
    for x in v.split("|"):
        f = (x.split(":") + ["", "", "", ""])[:4]
        out.append(f"{f[0]}:{f[1]} {f[2]}>{f[3]}")
    return sep.join(out)


def clinvar_combined_note(info, gt=""):
    """Step 6's note on a hit that is a combined allele of two neighbouring
    records, or part of one, from the tags bin/clinvar_join_adjacent.awk
    writes; bin/clinvar_hits.awk words it the same way. "." when none.
    On an unphased het call (0/1) the genotype cannot say the two parts share
    a haplotype, so the note says they may not."""
    a = gt.split("/")
    phase = ("; unphased, so the parts may sit on different haplotypes"
             if len(a) == 2 and all(x.isdigit() for x in a) and a[0] != a[1] else "")
    if info.get("PGP_JOINED_FROM"):
        return "combined allele of " + clinvar_variants(info["PGP_JOINED_FROM"], " + ") + phase
    if info.get("PGP_PART_OF"):
        return "part of a combined allele: " + clinvar_variants(info["PGP_PART_OF"], ", ") + phase
    if info.get("PGP_ONE_ALLELE"):
        return ("possibly one allele with " + clinvar_variants(info["PGP_ONE_ALLELE"], ", ")
                + ": check ClinVar for a combined record")
    return "."


def call_floor(path):
    """Step 6's call floor from its hits VCF header, CALL_FLOOR without the line."""
    floor = dict(CALL_FLOOR)
    try:
        with open_text(path) as f:
            for line in f:
                if not line.startswith("##"):
                    break
                if line.startswith("##PGP_CALL_FLOOR="):
                    for k, v in re.findall(r"(GQ|DP)=(\d+)", line):
                        floor[k] = int(v)
    except (OSError, EOFError, gzip.BadGzipFile, zlib.error) as e:
        raise Unreadable(f"{os.path.basename(path)}: {e}") from e
    return floor


def low_confidence(dp, gq, floor):
    """True when the call's GQ or DP is known and below the floor."""
    return (gq is not None and gq < floor["GQ"]) or (dp is not None and dp < floor["DP"])


def half_call(gt):
    """True for a GT of one alt allele and one missing one: 1/., ./1, 1|."""
    return bool(re.fullmatch(r"[1-9][/|]\.|\.[/|][1-9]", gt or ""))


def clinvar_zygosity(gt):
    """A ClinVar hit's genotype as bin/clinvar_hits.awk words it. One alt
    allele and one missing one (1/., ./1) is the half of a multiallelic call
    step 06 split with the other alt written as missing (a 1/2 at a poly-T
    tract), or a half call: never het, which would say the other copy is ref."""
    if half_call(gt):
        return f"one copy, other copy not called ref ({gt})"
    return {"0/1": "het", "1/0": "het", "0|1": "het", "1|0": "het", "1/1": "hom", "1|1": "hom"}.get(gt, gt)


def clinvar_row(r, floor=None):
    info = info_map(r[7]) if len(r) > 7 else {}
    rev = info.get("CLNREVSTAT", "")
    genes = [g.split(":")[0] for g in info.get("GENEINFO", "").split("|") if g]
    gt = r[9].split(":")[0] if len(r) > 9 else ""
    zyg = clinvar_zygosity(gt)
    fmt = dict(zip(r[8].split(":"), r[9].split(":"))) if len(r) > 9 else {}
    dp, gq = (int(fmt[k]) if fmt.get(k, "").isdigit() else None for k in ("DP", "GQ"))
    return {
        "stars": STARS.get(rev, 0),
        "chrom": r[0], "pos": int(r[1]), "ref": r[3], "alt": r[4],
        "genotype": zyg,
        "gene": ",".join(genes) or ".",
        "significance": info.get("CLNSIG", "").replace("_", " ") or ".",
        "review_status": rev.replace("_", " ") or ".",
        "combined_allele": clinvar_combined_note(info, gt),
        "dp": dp, "gq": gq,
        "low_confidence": low_confidence(dp, gq, floor or CALL_FLOOR),
    }


def clinvar_cdna_checks(d, s, p=None):
    """Step 06's cDNA read check (bin/clip_check.py), keyed by (chrom, pos,
    ref, alt), or the table at p. None when the file is absent: the Nextflow
    pipeline does not write it."""
    p = p or first_existing(d, [f"clinvar/{s}_clinvar_cdna_check.tsv"])
    if not p:
        return None
    try:
        rows = read_tsv(p)
    except (OSError, EOFError, gzip.BadGzipFile, zlib.error) as e:
        raise Unreadable(f"{os.path.basename(p)}: {e}") from e

    def num(v):
        return int(v) if (v or "").isdigit() else None
    out = {}
    for r in rows:
        if not (r.get("pos") or "").isdigit():
            continue
        out[(r.get("chrom"), int(r["pos"]), r.get("ref"), r.get("alt"))] = {
            "check": r.get("cdna_check") or "not_checked", "note": r.get("cdna_note") or "",
            "alt_reads": num(r.get("alt_reads")), "alt_clipped": num(r.get("alt_clipped")),
            "ref_reads": num(r.get("ref_reads")), "ref_clipped": num(r.get("ref_clipped")),
            "construct": r.get("construct_site") or None}
    return out


def clinvar_cdna_scan(d, s):
    """Step 06's genome-wide cDNA scan (bin/cdna_scan.py): {"note", "loci"},
    or None when the table is absent (a run before the scan existed)."""
    p = first_existing(d, [f"clinvar/{s}_clinvar_cdna_scan.tsv"])
    if not p:
        return None
    try:
        with open_text(p) as f:
            lines = f.read().splitlines()
    except (OSError, EOFError, gzip.BadGzipFile, zlib.error) as e:
        raise Unreadable(f"{os.path.basename(p)}: {e}") from e
    note = lines[0][len("# cdna_scan: "):] if lines and lines[0].startswith("# cdna_scan: ") else ""
    body = [ln for ln in lines if not ln.startswith("#")]
    loci = []
    for r in csv.DictReader(body, delimiter="\t"):
        if not (r.get("start") or "").isdigit():
            continue
        loci.append({"chrom": r.get("chrom"), "start": int(r["start"]), "end": int(r.get("end") or 0),
                     "junctions": int(r.get("junctions") or 0), "reads": int(r.get("reads") or 0),
                     "junction_list": r.get("junction_list") or ""})
    return {"note": note, "loci": loci}


def clinvar_deep_intronic(d, s):
    """Step 06's always-report list (the known deep intronic sites): one row
    per site with the sample's call, or None when the table is absent."""
    p = first_existing(d, [f"clinvar/{s}_clinvar_deep_intronic.tsv"])
    if not p:
        return None
    try:
        rows = read_tsv(p)
    except (OSError, EOFError, gzip.BadGzipFile, zlib.error) as e:
        raise Unreadable(f"{os.path.basename(p)}: {e}") from e
    return [{"chrom": r.get("chrom"), "pos": int(r["pos"]), "ref": r.get("ref"), "alt": r.get("alt"),
             "gene": r.get("gene") or ".", "hgvs": r.get("hgvs") or "", "clinvar_id": r.get("clinvar_id") or "",
             "call": r.get("call") or "no call", "detail": r.get("detail") or ""}
            for r in rows if (r.get("pos") or "").isdigit()]


def tabix_linear_index(path):
    """{contig: linear index} of a bgzipped VCF's .tbi: for each 16 kb window,
    the virtual file offset of the first record that can overlap it (standard
    library only, as the rest of this script)."""
    with gzip.open(path + ".tbi", "rb") as f:
        b = f.read()
    if b[:4] != b"TBI\x01":
        raise ValueError(f"{os.path.basename(path)}.tbi is not a tabix index")
    n_ref = struct.unpack_from("<i", b, 4)[0]
    l_nm = struct.unpack_from("<i", b, 32)[0]
    names = [x.decode() for x in b[36:36 + l_nm].split(b"\0")[:n_ref]]
    off, out = 36 + l_nm, {}
    for name in names:
        n_bin = struct.unpack_from("<i", b, off)[0]
        off += 4
        for _ in range(n_bin):
            n_chunk = struct.unpack_from("<i", b, off + 4)[0]
            off += 8 + 16 * n_chunk
        n_intv = struct.unpack_from("<i", b, off)[0]
        off += 4
        out[name] = struct.unpack_from(f"<{n_intv}Q", b, off)
        off += 8 * n_intv
    return out


def tabix_lines(path, index, chrom, pos):
    """The data lines of a bgzipped VCF that start at chrom:pos, read from the
    offset its tabix linear index gives for pos's window. A window with no
    offset recorded (0, or all ones) reads from the start of the file and
    skips the other contigs, which is slower but still right."""
    ioff = index.get(chrom)
    if not ioff:
        return []
    voff = ioff[min((pos - 1) >> 14, len(ioff) - 1)]
    if voff == 0xFFFFFFFFFFFFFFFF:
        voff = 0
    out, seen = [], False
    with open(path, "rb") as raw:
        raw.seek(voff >> 16)
        g = gzip.GzipFile(fileobj=raw)
        g.read(voff & 0xFFFF)
        for line in g:
            if line.startswith(b"#"):
                continue
            f = line.split(b"\t", 2)
            if len(f) < 3 or not f[1].isdigit():
                continue
            if f[0].decode() != chrom:
                if seen:
                    break
                continue
            seen = True
            if int(f[1]) > pos:
                break
            if int(f[1]) == pos:
                out.append(line.decode(errors="replace").rstrip("\n").split("\t"))
    return out


def vep_gnomad_af(vep, hits):
    """gnomAD exome and genome AF and the rsIDs of each hit from step 13's VEP
    output, as {(chrom, pos, ref, alt): (exome, genome, rsids)}; a hit with
    no record is left out. The fields are VEP's gnomADe_AF, gnomADg_AF and
    Existing_variation (--everything); rsids keeps the rs ids of
    Existing_variation ("rs123&COSV1" gives ["rs123"]), sorted, [] without
    the field."""
    fmt = None
    with open_text(vep) as f:
        for line in f:
            if not line.startswith("#"):
                break
            if line.startswith("##INFO=<ID=CSQ"):
                m = re.search(r"Format: ([^\"]+)", line)
                fmt = m.group(1).split("|") if m else None
    if not fmt or "Allele" not in fmt:
        return {}
    cols = [fmt.index(c) if c in fmt else None for c in ("gnomADe_AF", "gnomADg_AF")]
    ev = fmt.index("Existing_variation") if "Existing_variation" in fmt else None

    def num(v):
        vals = []
        for x in (v or "").split("&"):
            try:
                vals.append(float(x))
            except ValueError:
                pass
        return max(vals) if vals else None
    out = {}
    keys = sorted({(h["chrom"], h["pos"], h["ref"], h["alt"]) for h in hits})
    index = tabix_linear_index(vep) if keys else {}
    for key in keys:
        chrom, pos, ref, alt = key
        for r in tabix_lines(vep, index, chrom, pos):
            alts = r[4].split(",") if len(r) > 7 else []
            if r[3] != ref or alt not in alts:
                continue
            # VEP's Allele drops the base an indel's alleles share.
            snv = len(ref) == 1 and all(len(a) == 1 for a in alts)
            allele = alt if snv or any(a[:1] != ref[:1] for a in alts) else (alt[1:] or "-")
            got, ids = [None, None], set()
            for tr in info_map(r[7]).get("CSQ", "").split(","):
                v = tr.split("|")
                if len(v) != len(fmt) or (len(alts) > 1 and v[fmt.index("Allele")] != allele):
                    continue
                for i, c in enumerate(cols):
                    if c is not None and got[i] is None:
                        got[i] = num(v[c])
                if ev is not None:
                    ids.update(x.lower() for x in v[ev].split("&") if re.fullmatch(r"(?i)rs\d+", x))
            out[key] = (got[0], got[1], sorted(ids, key=lambda x: int(x[2:])))
    return out


def clinvar_frequencies(d, s, rows):
    """Put the gnomAD AF and the rsIDs (rsids, [] when unknown) from step
    13's VEP output on each ClinVar hit, and say where they came from or why
    there are none. The hits VCF's ID column is ClinVar's own variation id,
    not an rsID, so VEP is where a hit's rsID comes from."""
    vep = first_existing(d, [f"vep/{s}_vep.vcf.gz"])
    note = None
    if not vep:
        note = "no VEP output (step 13) in this folder"
    elif not os.path.exists(vep + ".tbi"):
        note = f"{os.path.relpath(vep, d)} has no .tbi index"
    found = {}
    if note is None:
        try:
            found = vep_gnomad_af(vep, rows)
        except (OSError, EOFError, ValueError, struct.error, gzip.BadGzipFile, zlib.error) as e:
            note = f"{os.path.relpath(vep, d)} could not be read: {e}"
    for x in rows:
        e, g, ids = found.get((x["chrom"], x["pos"], x["ref"], x["alt"]), (None, None, []))
        known = [v for v in (e, g) if v is not None]
        x["gnomad_exome_af"], x["gnomad_genome_af"] = e, g
        x["rsids"] = ids
        x["af"] = max(known) if known else None
        x["common"] = x["af"] is not None and x["af"] > COMMON_AF
    return {"source": os.path.relpath(vep, d) if note is None else None, "note": note, "common_af": COMMON_AF}


def sec_clinvar(d, s):
    p = first_existing(d, [f"clinvar/{s}_clinvar_hits.vcf"])
    if not p:
        return None, {}
    floor = call_floor(p)
    rows = [clinvar_row(r, floor) for r in vcf_records(p) if len(r) >= 8]
    rows.sort(key=lambda x: (-x["stars"], x["chrom"], x["pos"]))
    # The headline count and the stars leave the low-confidence calls out;
    # they stay in the hits, marked, and are counted on their own.
    sure = [x for x in rows if not x["low_confidence"]]
    by = {str(k): sum(1 for x in sure if x["stars"] == k) for k in (4, 3, 2, 1, 0)}
    data = {"count": len(sure), "by_stars": by, "hits": rows, "low_confidence": len(rows) - len(sure),
            "call_floor": floor, "frequency": clinvar_frequencies(d, s, rows)}
    # The cDNA read check, on each hit when step 06 wrote it.
    cdna = clinvar_cdna_checks(d, s)
    if cdna is not None:
        for x in rows:
            x["cdna"] = cdna.get((x["chrom"], x["pos"], x["ref"], x["alt"])) or {
                "check": "not_checked", "note": "not in the read-check table", "alt_reads": None,
                "alt_clipped": None, "ref_reads": None, "ref_clipped": None}
        data["cdna_flagged"] = sum(1 for x in rows if x["cdna"]["check"] == "flagged")
        data["construct_sites"] = sum(1 for x in rows if x["cdna"].get("construct"))
    scan = clinvar_cdna_scan(d, s)
    if scan is not None:
        data["cdna_scan"] = scan
    deep = clinvar_deep_intronic(d, s)
    if deep is not None:
        data["deep_intronic"] = deep
    return p, data


def pharmcat_json(d, s):
    return newest(d, [f"pharmcat/{s}.report.json", "pharmcat/*.report.json",
                      f"vcf/{s}.report.json", "vcf/*.report.json"])


def sec_pharmcat(d, s):
    p = pharmcat_json(d, s)
    if not p:
        return None, {}
    try:
        with open(p) as f:
            data = json.load(f)
    except (OSError, ValueError) as e:
        raise Unreadable(f"{os.path.basename(p)}: {e}")
    html = os.path.isfile(p[: -len(".json")] + ".html")
    return p, {"version": str(data.get("pharmcatVersion", "")) if isinstance(data, dict) else "",
               "html_report": html}


def sec_cpic(d, s):
    p = first_existing(d, [f"cpic/{s}_phenotypes.tsv"])
    if not p:
        return None, {}
    genes = []
    for row in read_tsv(p):
        genes.append({"gene": row.get("Gene", ""), "diplotype": row.get("Diplotype", ""),
                      "phenotype": row.get("Phenotype", ""), "status": row.get("Status", "")})
    rec = os.path.join(os.path.dirname(p), f"{s}_cpic_recommendations.txt")
    warnings, failed = [], False
    if os.path.isfile(rec):
        with open(rec, errors="replace") as f:
            for line in f:
                if "PARSING FAILED" in line:
                    failed = True
                if line.strip().startswith("WARNING:"):
                    warnings.append(line.strip()[len("WARNING:"):].strip())
    return p, {"genes": genes,
               "non_normal": sum(1 for g in genes if g["status"] == "non-normal"),
               "ambiguous": sum(1 for g in genes if g["status"] == "ambiguous"),
               "not_called": sum(1 for g in genes if g["status"] == "not called"),
               "unclassified": sum(1 for g in genes if g["status"] == "unclassified"),
               "parse_failed": failed or not genes,
               "warnings": warnings}


# The rows of the CYP2D6 depth check (bin/cyp2d6_depth_check.py) the reports
# show: its verdict and the MAPQ shares beside it.
DEPTH_CHECK_KEYS = ("status", "gene_mapq1_share", "gene_mapq20_share", "body_mapq1_share", "body_mapq20_share")


def sec_pypgx(d, s):
    p = first_existing(d, [f"pypgx/{s}_pypgx_summary.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    # A no-call by step 36's list (Indeterminate, FAILED, N/A...), in any case
    called = [r for r in rows if (r.get("Diplotype") or "").strip().lower() not in pgx_outside_calls.NO_CALL]
    cyp_row = next((r for r in rows if r.get("Gene") == "CYP2D6"), None)
    cyp = (cyp_row or {}).get("Diplotype", "")
    out = {"genes_total": len(rows), "genes_called": len(called), "cyp2d6": cyp or "not called"}
    # pypgx's copy-number call from read depth, one entry per BAM-based gene
    # (a VCF gene has none). A gene not callable on this reference gets its
    # own line below instead; one without a call ('N/A': pypgx failed, or
    # CYP2D6 after a failed depth check) says so, so no gene drops out.
    out["cnv_call"] = {}
    for r in rows:
        status = (r.get("Status") or "").strip()
        if (r.get("Source") or "").strip().upper() != "BAM" or status == pgx_outside_calls.NOT_CALLABLE:
            continue
        v = (r.get("CNV_call") or "").strip()
        if v not in ("", ".") and not v.startswith("N/A"):
            out["cnv_call"][r["Gene"]] = v
        elif status.lower() == "failed" or (r.get("Diplotype") or "").strip().upper() == "FAILED":
            out["cnv_call"][r["Gene"]] = "not called: pypgx failed"
        elif "depth check" in status or "depth check" in v:
            out["cnv_call"][r["Gene"]] = "not called: the CYP2D6 depth check failed"
        else:
            out["cnv_call"][r["Gene"]] = "not called: pypgx gave no copy-number call"
    out["not_callable"] = [r["Gene"] for r in rows
                           if (r.get("Status") or "").strip() == pgx_outside_calls.NOT_CALLABLE]
    out["low_gq"] = [{"gene": r["Gene"], "diplotype": (r.get("Diplotype") or "").strip(),
                      "min_gq": int(r["Min_GQ"]) if (r.get("Min_GQ") or "").strip().isdigit() else None,
                      "weakest": (r.get("Weakest_variant") or "").strip()}
                     for r in rows if (r.get("Quality") or "").strip() == pgx_outside_calls.LOW_GQ]
    if cyp_row is not None:
        out["cyp2d6_evidence"] = {"phenotype": (cyp_row.get("Phenotype") or "").strip(),
                                  "defining_filter": (cyp_row.get("Defining_filter") or "").strip(),
                                  "quality": (cyp_row.get("Quality") or "").strip()}
    dc = os.path.join(os.path.dirname(p), f"{s}_cyp2d6_depth_check.tsv")
    if os.path.isfile(dc):
        with open(dc, errors="replace") as f:
            kv = dict(line.rstrip("\n").split("\t", 1) for line in f if "\t" in line)
        out["depth_check"] = {k: kv[k] for k in DEPTH_CHECK_KEYS if k in kv}
    comp = os.path.join(os.path.dirname(p), f"{s}_pharmcat_comparison.tsv")
    if os.path.isfile(comp):
        crow = read_tsv(comp)
        out["comparison"] = {
            "conflicts": sum(1 for r in crow if r.get("Match") == "No"),
            "one_tool_only": sum(1 for r in crow if r.get("Match") in ("pypgx only", "PharmCAT only")),
            "concordant": sum(1 for r in crow if r.get("Match") == "Yes"),
        }
    return p, out


def sec_cyrius(d, s):
    p = first_existing(d, [f"cyrius/{s}_cyp2d6.tsv", "cyrius/*_cyp2d6.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    r = rows[0] if rows else {}
    return p, {"genotype": r.get("Genotype", "") or "none", "filter": r.get("Filter", "")}


def paralog_regions():
    """The Paraphase regions the report shows gene by gene, in the order of
    assets/paraphase_genes.tsv; the others are counted."""
    out = []
    for r in asset_rows("paraphase_genes.tsv"):
        if r["region"] not in out:
            out.append(r["region"])
    return out


def sec_paralogs(d, s):
    """Step 35b (Paraphase, PacBio HiFi): the genes of the regions of
    assets/paraphase_genes.tsv, a count of the other regions, the variant rows
    per region, and SMN1 beside step 35's Parascopy call when that table is
    there too."""
    p = first_existing(d, [f"paralogs/{s}_paraphase_copy_number.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    if rows and not {"region", "gene", "copy_number", "status"} <= set(rows[0]):
        raise Unreadable(f"{os.path.basename(p)} lacks the region, gene, copy_number or status column")
    PARALOG_REGIONS = paralog_regions()
    genes = [{"region": r["region"], "gene": r["gene"], "copy_number": r["copy_number"],
              "status": r["status"], "detail": r.get("detail") or "."}
             for r in rows if r["region"] in PARALOG_REGIONS]
    genes.sort(key=lambda g: PARALOG_REGIONS.index(g["region"]))
    others = [r for r in rows if r["region"] not in PARALOG_REGIONS]
    variants = {}
    vp = os.path.join(os.path.dirname(p), f"{s}_paraphase_variants.tsv")
    if os.path.isfile(vp):
        for r in read_tsv(vp):
            variants[r.get("region", "")] = variants.get(r.get("region", ""), 0) + 1
    smn1 = next((g["copy_number"] for g in genes if g["gene"] == "SMN1"), None)
    smn1 = int(smn1) if smn1 and smn1.isdigit() else None
    parascopy = None
    pp = os.path.join(os.path.dirname(p), f"{s}_smn_copy_number.tsv")
    if os.path.isfile(pp):
        parascopy = paraphase_report.parascopy_smn1(pp)
    return p, {"genes": genes,
               "other_regions": len(others),
               "other_regions_called": sum(1 for r in others if r["status"] == "called"),
               "variants": variants,
               "smn1": {"paraphase": smn1, "parascopy": parascopy,
                        "agree": paraphase_report.compare_smn1(smn1, parascopy)}}


# SMN2 exon 7, c.840 (the base that tells SMN2 from SMN1), GRCh38 1-based;
# SMN1's is paraphase_report.SMN1_EXON7.
SMN2_EXON7 = ("chr5", 70076526)


def parascopy_cn_at(rows, chrom, pos):
    """The copy number at one position from step 35's table: the
    paralog-specific value (psCN) of the row, or of its homologous region,
    that holds the position, with the row's aggregate value (agCN). As
    paraphase_report.parascopy_smn1: cn is None unless psCN's filter is PASS
    and its quality at least 20 (Parascopy's threshold)."""
    for r in rows:
        try:
            regions = [(r["chrom"], int(r["start"]) + 1, int(r["end"]))]
        except (KeyError, ValueError, TypeError):  # TypeError: a short row (None fields)
            continue
        for h in (r.get("homologous_regions") or "*").split(","):
            m = re.match(r"^([^:]+):(\d+)-(\d+)", h.strip())
            if m:
                regions.append((m.group(1), int(m.group(2)), int(m.group(3))))
        idx = next((i for i, (c, b, e) in enumerate(regions) if c == chrom and b <= pos <= e), None)
        if idx is None:
            continue
        cns = (r.get("psCN") or "").split(",")
        quals = (r.get("psCN_qual") or "").split(",")
        cn_txt = cns[idx] if idx < len(cns) else "?"
        q_txt = quals[idx] if idx < len(quals) else "*"
        filt = r.get("psCN_filter") or "."
        try:
            q = float(q_txt)
        except ValueError:
            q = None
        ok = filt == "PASS" and q is not None and q >= 20 and cn_txt.isdigit()
        return {"cn": int(cn_txt) if ok else None,
                "text": f"{cn_txt} (quality {q_txt}, filter {filt})" + ("" if ok else ", not resolved"),
                "agcn": r.get("agCN") or ".", "agcn_qual": r.get("agCN_qual") or ".",
                "agcn_filter": r.get("agCN_filter") or "."}
    return {"cn": None, "text": "no row covers exon 7", "agcn": ".", "agcn_qual": ".", "agcn_filter": "."}


def smn_reading(smn1):
    """What an SMN1 copy number means for SMA, in one line."""
    if smn1 is None:
        return "SMN1 not resolved (quality under 20 or filter not PASS): no carrier reading"
    if smn1 == 0:
        return "no SMN1 copy: the pattern of spinal muscular atrophy (SMA) itself"
    if smn1 == 1:
        return "one SMN1 copy: the pattern of an SMA carrier"
    if smn1 == 2:
        return ("two SMN1 copies: not a carrier by copy number; a 2+0 carrier (both copies on one chromosome) "
                "cannot be ruled out")
    return f"{smn1} SMN1 copies: not a carrier by copy number"


# Tables in paralogs/ that other parts of steps 35 and 35b write, not catalog rows.
NOT_CATALOG_TABLES = ("hba", "paraphase")
# agCN_filter and agCN_qual too: without them a catalog locus would read as
# "not resolved" instead of unreadable.
PARASCOPY_COLUMNS = {"chrom", "start", "end", "agCN", "agCN_filter", "agCN_qual", "psCN", "psCN_qual", "psCN_filter"}


def tsv_columns(path):
    """The header of a tab-separated table as a set ({} for an empty file), so
    a table with a header and no row is checked for its columns too."""
    with open_text(path) as f:
        return set(next(csv.reader(f, delimiter="\t"), []))


def parascopy_tables(d, s):
    """{table: path} of step 35's catalog tables (<s>_<table>_copy_number.tsv)."""
    out = {}
    for p in sorted(glob.glob(os.path.join(glob.escape(os.path.join(d, "paralogs", s)) + "_*_copy_number.tsv"))):
        name = os.path.basename(p)[len(s) + 1:-len("_copy_number.tsv")]
        if name not in NOT_CATALOG_TABLES and os.path.isfile(p):
            out[name] = p
    return out


def parascopy_locus(name, path, row):
    """One catalog table that is not SMN1's: the locus's aggregate copy number
    (agCN) per region of its profile, the catalog's genes, tag and reading.
    psCN is not read here: the catalog marks it reliable only for SMN1."""
    out = {"table": name, "genes": (row or {}).get("genes") or name, "tag": (row or {}).get("tag") or ".",
           "reading": (row or {}).get("reading") or "."}
    rows = read_tsv(path)
    if not PARASCOPY_COLUMNS <= tsv_columns(path):
        return dict(out, status="unreadable", cn=None,
                    text=f"unreadable: {os.path.basename(path)} lacks the copy-number columns")
    na = next((r["agCN_filter"] for r in rows if (r.get("agCN_filter") or "").startswith("not assessed")), None)
    if na:
        return dict(out, status="not assessed", cn=None, text=na)
    regions = []
    for r in rows:
        try:
            q = float(r.get("agCN_qual") or "")
        except ValueError:
            q = None
        cn = r.get("agCN") or "."
        ok = cn.isdigit() and r.get("agCN_filter") == "PASS" and q is not None and q >= 20
        regions.append((r.get("chrom"), r.get("start"), r.get("end"), cn, r.get("agCN_qual") or ".",
                        r.get("agCN_filter") or ".", ok))
    if not regions:
        return dict(out, status="not assessed", cn=None, text="not assessed: the table has no region row")
    values = {x[3] for x in regions if x[6]}
    if len(values) == 1 and all(x[6] for x in regions):
        cn = int(values.pop())
        quals = [float(x[4]) for x in regions]
        text = f"{cn} ({len(regions)} region(s), quality {min(quals):g} or more, filter PASS)"
        return dict(out, status="called", cn=cn, text=text)
    # every region, so a LowQual or discordant one is never hidden; no single
    # number, so no catalog reading either (the reading follows a call)
    text = "; ".join(f"{c}:{int(b) + 1 if str(b).isdigit() else b}-{e} {cn} (quality {q}, filter {f})"
                     for c, b, e, cn, q, f, _ in regions)
    called = any(x[6] for x in regions)
    return dict(out, status="called" if called else "not resolved", cn=None, reading=".",
                text=text + ("" if called else ": no region with a whole copy number, filter PASS and quality 20 or more"))


NO_SMN_TABLE = "no SMN1 table (step 35 wrote no SMN copy-number table)"


def sec_parascopy(d, s):
    """Step 35 (Parascopy, short reads): SMN1 and SMN2 copy number at exon 7,
    the aggregate copy number of the pair there and the SMA reading; the HBA
    count from depth; and every other catalog table, one locus each. Each part
    is read on its own, so a catalog without the SMN1 row still shows HBA and
    the other loci."""
    tables = parascopy_tables(d, s)
    hp = os.path.join(d, "paralogs", f"{s}_hba_copy_number.tsv")
    ip = os.path.join(d, "paralogs", f"{s}_cn_intervals.tsv")
    cp = os.path.join(d, "paralogs", f"{s}_cyp21a2_psv.tsv")
    dp = os.path.join(d, "paralogs", f"{s}_cn_dosage.tsv")
    p = (tables.get("smn") or next(iter(tables.values()), None)
         or next((x for x in (hp, ip, cp, dp) if os.path.isfile(x)), None))
    if not p:
        return None, {}
    if "smn" in tables:
        rows = read_tsv(tables["smn"])
        if not PARASCOPY_COLUMNS <= tsv_columns(tables["smn"]):
            raise Unreadable(f"{os.path.basename(tables['smn'])} lacks the chrom, start, end, agCN, agCN_filter, "
                             "agCN_qual, psCN, psCN_qual or psCN_filter column")
        smn1 = parascopy_cn_at(rows, *paraphase_report.SMN1_EXON7)
        smn2 = parascopy_cn_at(rows, *SMN2_EXON7)
        # Step 35 writes one row whose agCN_filter is "not assessed: <reason>"
        # when the locus was not run (no read over it, no model, Parascopy failed).
        na = next((r["agCN_filter"] for r in rows if (r.get("agCN_filter") or "").startswith("not assessed")), None)
        reading = (f"SMN1 and SMN2 {na}: no carrier reading" if na
                   else "no row of step 35's table covers SMN1 exon 7: no carrier reading"
                   if smn1["text"] == "no row covers exon 7" else smn_reading(smn1["cn"]))
    else:
        rows = []
        smn1 = {"cn": None, "text": NO_SMN_TABLE, "agcn": ".", "agcn_qual": ".", "agcn_filter": "."}
        smn2 = {"cn": None, "text": NO_SMN_TABLE}
        reading = f"{NO_SMN_TABLE}: no carrier reading"
    others = {n: t for n, t in tables.items() if n != "smn"}
    catalog = {r["table"]: r for r in asset_rows("paralog_loci.tsv")} if others else {}
    return p, {"smn1": smn1, "smn2": {"cn": smn2["cn"], "text": smn2["text"]},
               "agcn": smn1["agcn"], "agcn_qual": smn1["agcn_qual"], "agcn_filter": smn1["agcn_filter"],
               "rows": len(rows),
               # HBA1/HBA2 from depth (bin/hba_copy_number.py); None when
               # step 35 wrote no HBA table (a run from before the check)
               "hba": read_hba(hp) if os.path.isfile(hp) else None,
               # the interval rows of bin/cn_loci.tsv (recurrent CNV regions,
               # carrier deletions, AZF); None without step 35's table
               "cn_intervals": read_cn_intervals(ip) if os.path.isfile(ip) else None,
               # CYP21A2 from depth and PSV sites (bin/copy_number.py, the
               # CYP21A2 row): copies, module total, finding; None without it
               "cyp21a2": read_hba(cp) if os.path.isfile(cp) else None,
               # the gene-dosage rows (HP, LCE3BC, CCL3L1): total copies from
               # depth, one entry per locus; None without step 35's table
               "dosage": read_cn_dosage(dp) if os.path.isfile(dp) else None,
               "reading": reading,
               # every other catalog table, in the catalog's order (a table the
               # catalog does not name, from a PARALOG_LOCI of the user's, last)
               "loci": [parascopy_locus(n, others[n], catalog.get(n))
                        for n in sorted(others, key=lambda n: (list(catalog).index(n) if n in catalog
                                                               else len(catalog), n))]}


def read_cn_intervals(path):
    """Step 35's interval table (bin/copy_number.py intervals): the rows with
    a deletion or a duplication, and a count of the rows by status."""
    rows = read_tsv(path)
    cols = ("locus", "name", "chrom", "start", "end", "status", "copies", "expected", "change")
    if not set(cols) <= tsv_columns(path):
        raise Unreadable(f"{os.path.basename(path)} lacks the locus, name, status, copies, expected or change column")
    short_row(path, rows, cols)
    counts = {}
    for r in rows:
        counts[r["status"]] = counts.get(r["status"], 0) + 1
    return {"rows": len(rows), "by_status": counts,
            "changes": [{k: r[k] for k in ("locus", "name", "chrom", "start", "end", "copies", "expected", "change")}
                        for r in rows if r["change"] in ("deletion", "duplication")]}


DOSAGE_COLUMNS = ("locus", "status", "copies", "estimate", "reason", "reading")


def read_cn_dosage(path):
    """Step 35's gene-dosage table (HP, LCE3BC, CCL3L1): one entry per locus,
    in the table's order."""
    if not set(DOSAGE_COLUMNS) <= tsv_columns(path):
        raise Unreadable(f"{os.path.basename(path)} lacks the {', '.join(DOSAGE_COLUMNS[:-1])} or reading column")
    rows = read_tsv(path)
    short_row(path, rows, DOSAGE_COLUMNS)
    return [{k: r.get(k) or "." for k in DOSAGE_COLUMNS} for r in rows]


def short_row(path, rows, cols):
    """A row cut short of a column (csv gives None there) makes the table
    unreadable, so it never reads as a row of '.' values."""
    bad = next((i for i, r in enumerate(rows, 2) if any(r.get(k) is None for k in cols)), None)
    if bad:
        raise Unreadable(f"{os.path.basename(path)} line {bad} lacks a field")


def read_hba(path):
    """Step 35's HBA table (metric, value) as a dict, and its CYP21A2 table
    the same way; a table without a status row reads as not assessed, with
    that reason."""
    out = {}
    with open(path, errors="replace") as f:
        for line in f:
            c = line.rstrip("\n").split("\t", 1)
            if len(c) == 2 and c[0] != "metric":
                out[c[0]] = c[1]
    if out.get("status") not in ("called", "not assessed"):
        out = {"status": "not assessed", "message": f"not assessed: {os.path.basename(path)} has no status row"}
    return out


def sec_hla(d, s):
    p = first_existing(d, [f"hla_t1k/{s}_hla_genotype.tsv", f"hla/{s}_hla_genotype.tsv", "hla/*_hla_genotype.tsv"])
    if not p:
        return None, {}
    loci = []
    with open(p, errors="replace") as f:
        for line in f:
            c = line.rstrip("\n").split("\t")
            if len(c) < 3 or c[0].startswith("#"):
                continue
            # T1K: allele, abundance, quality, twice; '.', 0, -1 when there is no
            # second allele. Each quality stays with its allele.
            alleles, quals = [], []
            for i in (2, 5):
                a = c[i] if len(c) > i else ""
                if a and a != ".":
                    alleles.append(a)
                    quals.append(c[i + 2] if len(c) > i + 2 else "")
            # Step 36's rule (bin/pgx_outside_calls.py): an allele quality of 0 or
            # below withholds the whole gene; only HLA-A and HLA-B reach PharmCAT.
            low = any(hla_quality(q) <= 0 for q in quals)
            loci.append(hla_locus_extras({"gene": c[0], "alleles": alleles, "quality": quals, "low_confidence": low,
                                          "withheld_from_pharmcat": low and c[0] in pgx_outside_calls.HLA_GENES},
                                         c[1] if len(c) > 1 else ""))
    out = {"loci": loci, "derived": hla_derived(loci)}
    ot, kir = hla_optitype(d, s), kir_genotypes(d, s)
    if ot is not None:
        out["optitype"] = ot
    if kir is not None:
        out["kir"] = kir
    return p, out


def first_comment(path):
    """The first '#' line of a file step 08 wrote instead of a result, without the '#'."""
    with open(path, errors="replace") as f:
        for line in f:
            if line.startswith("#"):
                return line.lstrip("#").strip()
    return ""


def hla_optitype(d, s):
    """Step 08's second HLA caller (HLA_QC=true): per class I gene the row of
    its concordance table (pgx_outside_calls.py hla-concordance), or, when
    there is no table, why OptiType's result file says it did not run. None
    when OptiType was not asked for."""
    conc = first_existing(d, [f"hla_optitype/{s}_hla_concordance.tsv"])
    res = first_existing(d, [f"hla_optitype/{s}_optitype_result.tsv"])
    if not conc and not res:
        return None
    if not conc:
        return {"genes": [], "reason": first_comment(res) or "no concordance table: step 08 did not compare the callers"}
    genes = []
    with open(conc, errors="replace") as f:
        for r in csv.DictReader(f, delimiter="\t"):
            if (r.get("Gene") or "").strip():
                genes.append({"gene": r["Gene"].strip(), "t1k": (r.get("T1K") or "").strip(),
                              "optitype": (r.get("OptiType") or "").strip(), "agree": (r.get("Agree") or "").strip(),
                              "note": (r.get("Note") or "").strip()})
    return {"genes": genes, "reason": "" if genes else "the concordance table has no rows"}


def kir_genotypes(d, s):
    """Step 08's KIR pass (KIR=true): the genes T1K gave an allele, with each
    allele's quality, and the genes it gave none; or the line step 08 wrote
    for a BAM with too few KIR reads. None when KIR was not typed."""
    p = first_existing(d, [f"kir_t1k/{s}_kir_genotype.tsv", f"kir/{s}_kir_genotype.tsv"])
    if not p:
        return None
    typed, not_called = [], []
    with open(p, errors="replace") as f:
        for line in f:
            c = line.rstrip("\n").split("\t")
            if len(c) < 3 or c[0].startswith("#"):
                continue
            alleles, quals = [], []
            for i in (2, 5):
                if len(c) > i and c[i] not in ("", "."):
                    alleles.append(c[i])
                    quals.append(c[i + 2] if len(c) > i + 2 else "")
            if alleles:
                typed.append({"gene": c[0], "alleles": alleles, "quality": quals})
            else:
                not_called.append(c[0])
    reason = "" if typed or not_called else (first_comment(p) or "the KIR genotype file has no gene rows")
    return {"typed": typed, "not_called": not_called, "reason": reason}


def hla_quality(q):
    """T1K's quality as a number; one that is not a number counts as -1, as in step 36."""
    try:
        return float(q)
    except ValueError:
        return -1.0


# --- HLA: candidate caps, the HFE label, the DQ haplotypes and the allele table ---
# Rules in docs/08-hla-typing.md "What the report derives". The allele table is
# bin/hla_alleles.tsv (in bin/ because the report containers mount only bin/).

HLA_ALLELES_TSV = os.path.join(os.path.dirname(os.path.abspath(__file__)), "hla_alleles.tsv")
HLA_MAX_CANDIDATES = 3   # a quality-0 allele with more candidates is "no confident call"
HLA_HFE_LABEL = ("T1K sequence allele, not the HFE C282Y/H63D genotype: that genotype comes from the "
                 "variant calls (ClinVar section)")


def hla_name(x):
    """'HLA-DQB1*02:01:01' or 'DQB1*02:01:01' -> 'DQB1*02:01:01'; 'HLA-DQB1' -> 'DQB1'."""
    return x.strip()[4:] if x.strip().startswith("HLA-") else x.strip()


def hla_candidates(allele):
    """T1K's allele field: one allele, or equally likely candidates separated by commas."""
    return [a for a in (hla_name(x) for x in allele.split(",")) if a]


def hla_locus_extras(locus, n_diff):
    """Cap a quality-0 candidate list to HLA_MAX_CANDIDATES (the full list stays
    in the TSV), count the candidates, mark one homozygous call and label the
    HFE row."""
    counts, kept = [], []
    for a, q in zip(locus["alleles"], locus["quality"] + [""] * len(locus["alleles"])):
        cands = hla_candidates(a)
        counts.append(len(cands))
        kept.append(",".join(cands[:HLA_MAX_CANDIDATES]) if hla_quality(q) <= 0 and len(cands) > HLA_MAX_CANDIDATES
                    else a)
    locus["alleles"] = kept
    locus["candidates"] = counts
    # T1K writes one allele and '. 0 -1' for a homozygous gene (or a second allele it did not see).
    locus["one_allele_type"] = n_diff.strip() == "1" and len(kept) == 1
    if hla_name(locus["gene"]) == "HFE":
        locus["label"] = HLA_HFE_LABEL
    return locus


def hla_match(allele, group):
    """'yes' when every candidate of the allele is in the group (B*57:01 holds
    B*57:01:01:02), 'no' when none is, 'mixed' when only some are."""
    hits = [c == group or c.startswith(group + ":") for c in hla_candidates(allele)]
    return "yes" if hits and all(hits) else ("mixed" if any(hits) else "no")


def hla_read_table(path=HLA_ALLELES_TSV):
    """bin/hla_alleles.tsv as a list of dicts (tests/test_hla_report.py checks its schema)."""
    rows, head = [], None
    with open(path) as f:
        for line in f:
            if line.startswith("#") or not line.strip():
                continue
            c = line.rstrip("\n").split("\t")
            if head is None:
                head = c
                continue
            rows.append(dict(zip(head, c)))
    return rows


def hla_derived(loci, table=None):
    """The DQ haplotypes, the DQB1*02 dose and each table allele present or absent."""
    by = {}
    for l in loci:
        by.setdefault(hla_name(l["gene"]), l)

    def usable(gene):
        l = by.get(gene)
        if not l or not l["alleles"]:
            return None, f"HLA-{gene} not typed"
        if l["low_confidence"]:
            return None, f"HLA-{gene} low confidence"
        return l, ""

    def carries(l, group):
        """'yes', 'no' or 'mixed' over both alleles of a locus."""
        m = [hla_match(a, group) for a in l["alleles"]]
        return "yes" if "yes" in m else ("mixed" if "mixed" in m else "no")

    def both(a, b):
        return "present" if a == b == "yes" else ("absent" if "no" in (a, b) else "uncertain")

    dqa, ra = usable("DQA1")
    dqb, rb = usable("DQB1")
    coeliac = {"status": "not_determined", "reason": "; ".join(r for r in (ra, rb) if r),
               "dq2_5": None, "dq2_5_phase": None, "dq2_2": None, "dq8": None, "dqb1_02_dose": None}
    if dqb:
        m = [hla_match(a, "DQB1*02") for a in dqb["alleles"]]
        if "mixed" not in m:
            n = m.count("yes")
            coeliac["dqb1_02_dose"] = 2 * n if dqb["one_allele_type"] else n
    if dqa and dqb:
        coeliac.update(status="derived", reason="")
        coeliac["dq2_5"] = both(carries(dqa, "DQA1*05"), carries(dqb, "DQB1*02"))
        coeliac["dq2_2"] = both(carries(dqa, "DQA1*02:01"), carries(dqb, "DQB1*02:02"))
        coeliac["dq8"] = both(carries(dqa, "DQA1*03"), carries(dqb, "DQB1*03:02"))
        if coeliac["dq2_5"] == "present":
            drb1, why = usable("DRB1")
            dr3 = carries(drb1, "DRB1*03") if drb1 else None
            coeliac["dq2_5_phase"] = ("likely cis (DRB1*03 present)" if dr3 == "yes"
                                      else "possible trans (no DRB1*03)" if dr3 == "no"
                                      else f"phase unknown ({why or 'DRB1*03 uncertain'})")
    alleles = []
    for r in (hla_read_table() if table is None else table):
        gene = r["allele"].split("*")[0]
        l, why = usable(gene)
        state = carries(l, r["allele"]) if l else None
        alleles.append({"allele": r["allele"], "kind": r["kind"], "association": r["association"],
                        "pmids": r["pmids"], "citation": r["citation"],
                        "state": {"yes": "present", "no": "absent", "mixed": "uncertain", None: "not_determined"}[state],
                        "reason": why})
    return {"coeliac": coeliac, "alleles": alleles}


def sec_prs(d, s):
    p = first_existing(d, [f"prs/{s}_prs_summary.tsv"])
    if not p:
        return None, {}
    rows = []
    for r in read_tsv(p):
        pct = r.get("Percentile", "") or ""
        rows.append({"condition": r.get("Condition", ""), "pgs_id": r.get("PGS_ID", ""),
                     "score": r.get("Score_SUM", ""), "matched": r.get("Variants_Matched", ""),
                     "total": r.get("Variants_Total", ""),
                     "chrx": r.get("ChrX", ""),
                     # "12 dropped", "12 kept", "12 dropped, strand check failed", "0"; empty in older tables
                     "ambiguous": r.get("Ambiguous", ""),
                     # the sex a sex-specific score applies to, and the check against the declared sex
                     "applies_to": r.get("Applies_To", ""), "sex_check": r.get("Sex_Check", ""),
                     "percentile": "" if pct == "NA" else pct,
                     "group": "" if r.get("Ancestry_Group", "NA") == "NA" else r.get("Ancestry_Group", ""),
                     # the groups the PGS Catalog has an OR or HR of the score for, and the
                     # one of the sample's group ("HR 1.53 (EUR, PPM000591)"); empty in older tables
                     # None when the table has no Evaluated_In column (a score list without
                     # evaluations, or an older table): nothing is known, so nothing is said
                     "evaluated_in": ([x for x in r["Evaluated_In"].split(";") if x]
                                      if r.get("Evaluated_In") is not None else None),
                     "effect": r.get("Effect", "") or ""})
    data = {"scores": rows,
            # Percentiles exist only when step 25 ran pgsc_calc with an ancestry panel.
            "adjusted": any(r["percentile"] for r in rows)}
    a = first_existing(d, [f"ancestry/{s}_ancestry.tsv", f"prs/{s}_ancestry.tsv"])
    if a:
        kv = {r.get("key"): r.get("value", "") for r in read_tsv(a) if r.get("key")}
        data["ancestry"] = {"population": kv.get("population", ""), "panel": kv.get("reference_panel", ""),
                            "low_confidence": kv.get("population_low_confidence", "")}
    return p, data


# --- PRS with pgsc_calc (steps 25 and 26, and the PRS processes) ------------------
# The PGS Catalog's pgsc_calc scores the sample. Step 25 and the Nextflow
# PRS_PREPARE and PRS_SUMMARY processes run the two commands below, so the
# score files pgsc_calc reads and the table the reports read are made once.

# Input stays last: readers take it as the last column.
PRS_COLUMNS = ["Condition", "PGS_ID", "Score_SUM", "Variants_Matched", "Variants_Total",
               "Matched_Pct", "Percentile", "Ancestry_Group", "ChrX", "Ambiguous", "Applies_To",
               "Sex_Check", "Input"]
# Written before Input when the score list names PGS Catalog evaluations
# (ppm_id, ancestry and effect in assets/pgs_scores.tsv): the groups a score
# was evaluated in, and the effect measured in the sample's group.
PRS_EFFECT_COLUMNS = ["Evaluated_In", "Effect"]
# prs-format writes it beside the scores: per score, its chrX rows and
# whether they were kept (scored) or left out (no sex given).
CHRX_FILE = "chrx_rows.tsv"
# prs-format writes it beside the scores too: per score, the sex it applies
# to (applies_to of assets/pgs_scores.tsv) and its strand-ambiguous rows (an
# A/T or C/G pair) with what happened to them.
SCORE_NOTES_FILE = "score_notes.tsv"
# prs-format writes it beside the scores as well: per score, the PGS Catalog
# evaluations of assets/pgs_scores.tsv (pgs_id, ppm_id, ancestry, effect).
EVALS_FILE = "score_evaluations.tsv"
SEXES = ("female", "male")
# The 1000 Genomes groups pgsc_calc places a sample in, which the ancestry
# column of assets/pgs_scores.tsv names.
PRS_GROUPS = ("AFR", "AMR", "EAS", "EUR", "SAS")
# Sex_Check of a sex-specific score: the declared sex is the one it applies
# to, another one (no percentile then), or none was declared.
SEX_APPLIES = "applies"
SEX_NOT_APPLICABLE = "not applicable to the declared sex"
SEX_NOT_DECLARED = "no sex declared"
COMPLEMENT = {"A": "T", "T": "A", "C": "G", "G": "C"}
# Columns of a scoring file that make a score non-additive: pgsc_calc scores
# them, but this pipeline's custom file keeps effect_weight only.
NON_ADDITIVE = ("dosage_0_weight", "dosage_1_weight", "dosage_2_weight")
AUTOSOMES = frozenset(str(i) for i in range(1, 23))


def pgs_id_of(path):
    """PGS000018 for PGS000018.txt.gz, PGS000018_hmPOS_GRCh38.txt.gz and the like."""
    name = os.path.basename(path)
    for suffix in (".gz", ".txt", ".tsv"):
        if name.endswith(suffix):
            name = name[: -len(suffix)]
    return name.split("_hmPOS_")[0]


def read_labels(path):
    """ID -> trait label of assets/pgs_scores.tsv (columns pgs_id, trait_reported)."""
    out = {}
    if not path:
        return out
    with open(path, encoding="utf-8") as f:
        for r in csv.DictReader((line for line in f if not line.startswith("#")), delimiter="\t"):
            if r.get("pgs_id"):
                out[r["pgs_id"].strip()] = (r.get("trait_reported") or "").strip()
    return out


def read_applies_to(path):
    """ID -> the sex a score applies to, from the optional applies_to column of
    assets/pgs_scores.tsv: 'female', 'male', or '' for a score of both sexes.
    Raises Unreadable for any other value."""
    out = {}
    if not path:
        return out
    with open(path, encoding="utf-8") as f:
        for r in csv.DictReader((line for line in f if not line.startswith("#")), delimiter="\t"):
            pid = (r.get("pgs_id") or "").strip()
            sex = (r.get("applies_to") or "").strip().lower()
            if not pid:
                continue
            if sex and sex not in SEXES:
                raise Unreadable(f"{os.path.basename(path)}: applies_to of {pid} is '{sex}', expected female, male or empty")
            out[pid] = sex
    return out


def read_metrics(path):
    """ID -> [(ancestry, effect, ppm_id), ...]: the PGS Catalog evaluations of
    each score, from the optional ppm_id, ancestry and effect columns of
    assets/pgs_scores.tsv (entries separated by ";", in the same order in the
    three). Raises Unreadable when a row with evaluations has no pgs_id or
    repeats one, the three lists differ in length, an ancestry is not a 1000
    Genomes group or appears twice, or an effect is not "OR <number>" or
    "HR <number>"."""
    out = {}
    if not path:
        return out
    with open(path, encoding="utf-8") as f:
        for r in csv.DictReader((line for line in f if not line.startswith("#")), delimiter="\t"):
            pid = (r.get("pgs_id") or "").strip()
            cols = [[x.strip() for x in (r.get(k) or "").split(";")] if (r.get(k) or "").strip() else []
                    for k in ("ppm_id", "ancestry", "effect")]
            name = os.path.basename(path)
            if not any(cols):
                continue
            if not pid:
                raise Unreadable(f"{name}: a row has ppm_id, ancestry or effect entries but no pgs_id")
            if pid in out:
                raise Unreadable(f"{name}: {pid} has evaluations on two rows")
            ppm, anc, eff = cols
            if not (len(ppm) == len(anc) == len(eff)):
                raise Unreadable(f"{name}: {pid} has {len(ppm)} ppm_id, {len(anc)} ancestry and {len(eff)} effect entries")
            for p_, a_, e_ in zip(ppm, anc, eff):
                if a_ not in PRS_GROUPS:
                    raise Unreadable(f"{name}: ancestry of {pid} is '{a_}', expected one of {', '.join(PRS_GROUPS)}")
                if not re.fullmatch(r"PPM\d{6}", p_):
                    raise Unreadable(f"{name}: ppm_id of {pid} is '{p_}', expected PPM and six digits")
                if not re.fullmatch(r"(OR|HR) \d+(\.\d+)?", e_):
                    raise Unreadable(f"{name}: effect of {pid} is '{e_}', expected 'OR <number>' or 'HR <number>'")
            if len(set(anc)) != len(anc):
                raise Unreadable(f"{name}: {pid} lists an ancestry twice")
            out[pid] = list(zip(anc, eff, ppm))
    return out


class Fasta:
    """Bases of an uncompressed FASTA with its .fai, read by seeking."""

    def __init__(self, path):
        self.idx = {}
        try:
            with open(path + ".fai") as f:
                for line in f:
                    p = line.rstrip("\n").split("\t")
                    if len(p) >= 5:
                        self.idx[p[0]] = tuple(int(x) for x in p[1:5])
            self.f = open(path, "rb")
        except (OSError, ValueError) as e:
            raise Unreadable(f"cannot read the reference {path} and its .fai: {e}")
        if self.f.read(2) == b"\x1f\x8b":
            raise Unreadable(f"{path} is compressed; the strand check reads an uncompressed FASTA")

    def base(self, chrom, pos):
        """The base at 1-based POS of CHROM (with or without chr), upper case;
        '' when the contig is not in the reference or POS is outside it."""
        c = self.idx.get(chrom) or self.idx.get("chr" + chrom)
        if not c or not 1 <= pos <= c[0]:
            return ""
        length, offset, bases, width = c
        i = pos - 1
        self.f.seek(offset + (i // bases) * width + i % bases)
        return self.f.read(1).decode("ascii", "replace").upper()


def ambiguous_pair(ea, oa):
    """True for a strand-ambiguous SNV: its two alleles are complements (A/T, C/G)."""
    return len(ea) == 1 and len(oa) == 1 and COMPLEMENT.get(ea) == oa


def strand_check(rows, ref):
    """Whether a score's alleles are given on the forward strand, so its
    strand-ambiguous rows can be kept: (passed, rows that agree with the
    reference, rows that agree only after a strand flip). ROWS are (chrom,
    pos, effect, other) of the rows written. A row counts when both alleles
    are single bases and not an ambiguous pair: it agrees when the reference
    base is one of its alleles, and is flipped when the base is the
    complement of one of them instead. A score passes with at least one
    agreeing row and no flipped one; a row whose contig is not in the
    reference is not counted."""
    agree = flipped = 0
    for chrom, pos, ea, oa in rows:
        if len(ea) != 1 or len(oa) != 1 or ambiguous_pair(ea, oa) or not pos.isdigit():
            continue
        b = ref.base(chrom, int(pos))
        if b not in COMPLEMENT:
            continue
        if b in (ea, oa):
            agree += 1
        elif b in (COMPLEMENT.get(ea), COMPLEMENT.get(oa)):
            flipped += 1
    return agree > 0 and flipped == 0, agree, flipped


def single_allele(a):
    a = (a or "").strip().upper()
    return a if a and "/" not in a and a != "." else ""


def format_pgs(path, out_dir, label=None, keep_x=False, ref=None):
    """Write the GRCh38-harmonised scoring file PATH as a custom GRCh38 file
    pgsc_calc reads (chr_name and chr_position from hm_chr and hm_pos), and
    return (pgs_id, rows written, rows dropped off the autosomes,
    [(chrom, pos, allele), ...], chrX rows, ambiguous) with the effect and
    other alleles of every row written. Chromosomes 1 to 22 are kept, and chrX
    with keep_x (the sample's sex is known): pgsc_calc converts the sample
    with plink2, which stops on a chrX record when no sex is given. chrY and
    MT rows are always dropped. AMBIGUOUS is (strand-ambiguous rows, what
    happens to them, agreeing rows, flipped rows): without REF they are
    written and pgsc_calc drops them ("dropped"); with REF (keep_ambiguous)
    they are written and kept when the score passes strand_check ("kept"),
    and left out of the file when it does not ("dropped, strand check
    failed"). Raises Unreadable for a file that is not harmonised to GRCh38,
    lacks a column, is not additive or has no row."""
    header, cols, rows, alleles = {}, None, [], []
    with open_text(path) as f:
        for line in f:
            line = line.rstrip("\n")
            if line.startswith("#"):
                k, _, v = line[1:].partition("=")
                header[k.strip()] = v.strip()
                continue
            if cols is None:
                cols = line.split("\t")
                continue
            if line:
                rows.append(line.split("\t"))
    name = os.path.basename(path)
    if header.get("HmPOS_build") != "GRCh38":
        raise Unreadable(f"{name} has #HmPOS_build='{header.get('HmPOS_build', '')}', expected GRCh38 "
                         "(use the PGS Catalog's Harmonized/<id>_hmPOS_GRCh38 file)")
    c = {k: i for i, k in enumerate(cols or [])}
    missing = [k for k in ("hm_chr", "hm_pos", "effect_allele", "effect_weight") if k not in c]
    if missing:
        raise Unreadable(f"{name} has no {', '.join(missing)} column")
    if any(k in c for k in NON_ADDITIVE):
        raise Unreadable(f"{name} has dosage weights (a non-additive score); only additive scores are supported")

    def get(r, k):
        return r[c[k]].strip() if k in c and c[k] < len(r) else ""

    for k in ("is_dominant", "is_recessive"):
        if k in c and any(get(r, k) == "True" for r in rows):
            raise Unreadable(f"{name} has {k}=True rows (a non-additive score); only additive scores are supported")
    pid = header.get("pgs_id") or pgs_id_of(path)
    trait = (label or header.get("trait_reported") or pid).replace("=", "-").replace("\n", " ")
    out = os.path.join(out_dir, f"{pid}.txt.gz")
    counts = {"off": 0, "x": 0}

    def placed():
        """(chrom, pos, effect, other, weight) of each row pgsc_calc gets; one
        pass over ROWS each call, so no second copy of a large score is held."""
        counts["off"] = counts["x"] = 0
        for r in rows:
            chrom, pos, ea, ew = get(r, "hm_chr"), get(r, "hm_pos"), get(r, "effect_allele").upper(), get(r, "effect_weight")
            if not (chrom and pos and ea and ew):
                continue   # a row the catalog could not place on GRCh38
            chrom = chrom[3:] if chrom.startswith("chr") else chrom
            if chrom == "X":
                counts["x"] += 1
            if chrom not in AUTOSOMES and not (keep_x and chrom == "X"):
                counts["off"] += 1   # chrY, MT, and chrX unless the sex is known: plink2 refuses chrX without it
                continue
            yield chrom, pos, ea, single_allele(get(r, "other_allele")) or single_allele(get(r, "hm_inferOtherAllele")), ew

    agree = flipped = 0
    amb, drop_amb = "dropped", False
    if ref is not None:
        passed, agree, flipped = strand_check(((c, p, ea, oa) for c, p, ea, oa, _ in placed()), ref)
        if passed:
            amb = "kept"
        else:
            # pgsc_calc's keep_ambiguous applies to every score of the run,
            # so a score that fails loses its ambiguous rows here.
            amb, drop_amb = "dropped, strand check failed", True
    n = n_amb = 0
    with gzip.open(out + ".tmp", "wt", encoding="utf-8") as w:
        w.write(f"#pgs_id={pid}\n#pgs_name={pid}\n#trait_reported={trait}\n#genome_build=GRCh38\n")
        w.write("chr_name\tchr_position\teffect_allele\tother_allele\teffect_weight\n")
        for chrom, pos, ea, oa, ew in placed():
            if ambiguous_pair(ea, oa):
                n_amb += 1
                if drop_amb:
                    continue
            w.write(f"{chrom}\t{pos}\t{ea}\t{oa}\t{ew}\n")
            n += 1
            alleles.append((f"chr{chrom}", pos, ea))
            if oa:
                alleles.append((f"chr{chrom}", pos, oa))
    if n == 0:
        os.remove(out + ".tmp")
        raise Unreadable(f"{name} has no autosomal row with hm_chr, hm_pos, effect_allele and effect_weight")
    os.replace(out + ".tmp", out)
    return pid, n, counts["off"], alleles, counts["x"], (n_amb, amb, agree, flipped)


def prs_format_main(argv):
    ap = argparse.ArgumentParser(prog="collect_summary.py prs-format",
                                 description="Write the scoring files as the custom GRCh38 files pgsc_calc reads.")
    ap.add_argument("--scores", required=True, help="folder of GRCh38-harmonised PGS Catalog scoring files")
    ap.add_argument("--ids", default="", help="comma-separated PGS ids to use (default: every file in --scores)")
    ap.add_argument("--labels", help="assets/pgs_scores.tsv: the trait label of each id")
    ap.add_argument("--out", required=True, help="folder for the files pgsc_calc reads")
    ap.add_argument("--alleles", required=True,
                    help="written: CHROM POS ALLELE of every effect and other allele, sorted and unique")
    ap.add_argument("--keep-x", action="store_true",
                    help="keep the chrX rows (the sample's sex is known and reaches pgsc_calc's plink2); "
                         "without it they are dropped with chrY and MT")
    ap.add_argument("--keep-ambiguous", action="store_true",
                    help="pgsc_calc runs with --keep_ambiguous: keep each score's strand-ambiguous rows (A/T, C/G) "
                         "when its other rows agree with --reference on the forward strand, and leave them out of "
                         "a score that fails that check")
    ap.add_argument("--reference", help="the uncompressed GRCh38 FASTA with its .fai, for --keep-ambiguous")
    a = ap.parse_args(argv)
    if a.keep_ambiguous and not a.reference:
        print("ERROR: --keep-ambiguous needs --reference, the FASTA the strand check reads", file=sys.stderr)
        return 1
    try:
        applies = read_applies_to(a.labels)
        metrics = read_metrics(a.labels)
        ref = Fasta(a.reference) if a.keep_ambiguous else None
    except (Unreadable, OSError) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        return 1
    labels = read_labels(a.labels)
    files = sorted(glob.glob(os.path.join(a.scores, "*.txt.gz")) + glob.glob(os.path.join(a.scores, "*.txt")))
    if a.ids:
        want = [x.strip() for x in a.ids.split(",") if x.strip()]
        by_id = {pgs_id_of(p): p for p in files}
        missing = [x for x in want if x not in by_id]
        if missing:
            print(f"ERROR: no scoring file for {', '.join(missing)} in {a.scores}", file=sys.stderr)
            return 1
        files = [by_id[x] for x in want]
    if not files:
        print(f"ERROR: no scoring file (*.txt.gz or *.txt) in {a.scores}", file=sys.stderr)
        return 1
    os.makedirs(a.out, exist_ok=True)
    allele_set, chrx, notes = set(), [], []
    try:
        for p in files:
            pid, n, off, al, nx, (n_amb, amb, agree, flipped) = format_pgs(
                p, a.out, labels.get(pgs_id_of(p)), a.keep_x, ref)
            allele_set.update(al)
            chrx.append((pid, nx))
            notes.append((pid, applies.get(pid, ""), n_amb, amb))
            kept = f", {nx} on chrX scored" if a.keep_x and nx else ""
            dropped = f", {off} off the autosomes dropped" if off else ""
            ambig = f", {n_amb} strand-ambiguous {amb}" if n_amb else ""
            if ref is not None:
                ambig += f" (strand check: {agree} rows agree with the reference, {flipped} flipped)"
            sex = f", {applies[pid]} only" if applies.get(pid) else ""
            print(f"  {pid}: {n} GRCh38 rows{kept}{dropped}{ambig}{sex} ({labels.get(pid) or 'label from the file'})")
    except (Unreadable, OSError, EOFError, zlib.error) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        return 1
    with open(os.path.join(a.out, CHRX_FILE), "w") as w:
        w.write("pgs_id\tchrx_rows\tchrx\n")
        for pid, nx in chrx:
            w.write(f"{pid}\t{nx}\t{'scored' if a.keep_x else 'left out'}\n")
    with open(os.path.join(a.out, SCORE_NOTES_FILE), "w") as w:
        w.write("pgs_id\tapplies_to\tambiguous_rows\tambiguous\n")
        for pid, sex, n_amb, amb in notes:
            w.write(f"{pid}\t{sex}\t{n_amb}\t{amb}\n")
    with open(os.path.join(a.out, EVALS_FILE), "w", encoding="utf-8") as w:
        w.write("pgs_id\tppm_id\tancestry\teffect\n")
        for pid, *_ in notes:
            ev = metrics.get(pid)
            if ev:
                w.write("\t".join([pid] + [";".join(e[i] for e in ev) for i in (2, 0, 1)]) + "\n")

    def key(t):
        return (t[0], int(t[1]) if t[1].isdigit() else 0, t[2])
    with open(a.alleles + ".tmp", "w") as w:
        for t in sorted(allele_set, key=key):
            w.write("\t".join(t) + "\n")
    os.replace(a.alleles + ".tmp", a.alleles)
    return 0


def fmt_num(x):
    """A number as pgsc_calc wrote it, without a trailing .0: 62.0 -> 62."""
    try:
        return f"{float(x):.10g}"
    except (TypeError, ValueError):
        return "NA"


def read_gz_tsv(path):
    with open_text(path) as f:
        return list(csv.DictReader(f, delimiter="\t"))


# pgscatalog-match's log line for a score under the minimum overlap.
BELOW_RE = re.compile(r"Score (\S+) fails minimum matching threshold \(([0-9.]+)% variants match\)")


def chrx_matched(results, sampleset):
    """PGS id -> chrX rows pgsc_calc matched, from its match log; None when
    there is no log. In a run with an ancestry panel pgsc_calc keeps only
    the variants it found in the panel too (the rest are 'excluded'), so a
    chrX row the panel lacks is not scored even with the sex."""
    logs = sorted(glob.glob(os.path.join(results, sampleset, "match", "*_log.csv.gz")))
    if not logs:
        return None
    out = {}
    for p in logs:
        with open_text(p) as f:
            head = next(csv.reader([f.readline()]), [])
            c = {k: i for i, k in enumerate(head)}
            if not {"accession", "chr_name", "match_status"} <= set(c):
                return None
            for line in f:
                # Most lines are not on chrX; split only those that may be.
                if ",X," not in line and ",23," not in line:
                    continue
                r = next(csv.reader([line]))
                if len(r) > max(c.values()) and r[c["chr_name"]] in ("X", "23") and r[c["match_status"]] == "matched":
                    pid = pgs_id_of(r[c["accession"]])
                    out[pid] = out.get(pid, 0) + 1
    return out


def prs_table(results, sampleset, scores_dir, input_kind, zero_matches=False, below_log=None, sex="", panel=False):
    """Rows of the summary table (PRS_COLUMNS) and the ancestry rows (key,
    value), from pgsc_calc's output folder RESULTS for SAMPLESET. Totals and
    labels come from the files prs-format wrote in SCORES_DIR. With
    zero_matches (pgsc_calc stopped because no score variant is in the
    sample's genotypes) every score is reported unmatched. With below_log
    (pgsc_calc stopped because every score matched under its minimum
    overlap, so it published no match summary) no score has a sum, the
    matched count is unknown (NA), and the match rate is read from that log
    when pgscatalog-match printed it. SEX is the declared sex ('' for none):
    a score that applies to the other sex keeps its sum and loses its
    percentile. PANEL: pgsc_calc ran with an ancestry panel, so a chrX row
    counts as scored only when its match log matched it."""
    totals, labels = {}, {}
    for p in sorted(glob.glob(os.path.join(scores_dir, "*.txt.gz"))):
        pid, n = pgs_id_of(p), 0
        with open_text(p) as f:
            for line in f:
                if line.startswith("#trait_reported="):
                    labels[pid] = line.rstrip("\n").split("=", 1)[1]
                elif not line.startswith("#"):
                    n += 1
        totals[pid] = n - 1   # minus the column header
    if not totals:
        raise Unreadable(f"no formatted scoring file in {scores_dir}")
    # "8 scored", "8 left out" or "0"; empty for scores formatted before the file existed
    chrx, chrx_scored = {}, set()
    if os.path.isfile(os.path.join(scores_dir, CHRX_FILE)):
        for r in read_tsv(os.path.join(scores_dir, CHRX_FILE)):
            n = r.get("chrx_rows", "")
            chrx[r.get("pgs_id", "")] = f"{n} {r.get('chrx', '')}" if n not in ("", "0") else n
            if n not in ("", "0") and r.get("chrx") == "scored":
                chrx_scored.add(r.get("pgs_id", ""))
    # "" for a score of both sexes; empty for scores formatted before the file existed
    applies, ambig, evals = {}, {}, {}
    notes_path = os.path.join(scores_dir, SCORE_NOTES_FILE)
    if os.path.isfile(notes_path):
        for r in read_tsv(notes_path):
            pid, n = r.get("pgs_id", ""), r.get("ambiguous_rows", "")
            applies[pid] = r.get("applies_to", "")
            ambig[pid] = f"{n} {r.get('ambiguous', '')}" if n not in ("", "0") else n
    # The evaluations prs-format copied from assets/pgs_scores.tsv; none in a
    # folder formatted before it wrote them.
    if os.path.isfile(os.path.join(scores_dir, EVALS_FILE)):
        evals = read_metrics(os.path.join(scores_dir, EVALS_FILE))
    sex = (sex or "").lower()
    matched, sums, pct, group = {}, {}, {}, {}
    pops, rate = {}, {}
    if below_log:
        with open(below_log, errors="replace") as f:
            for m in BELOW_RE.finditer(f.read()):
                rate[pgs_id_of(m.group(1))] = f"{float(m.group(2)):.1f}"
    elif not zero_matches:
        summ = os.path.join(results, sampleset, "match", f"{sampleset}_summary.csv")
        if not os.path.isfile(summ):
            raise Unreadable(f"pgsc_calc wrote no match summary ({summ})")
        with open(summ, newline="") as f:
            for r in csv.DictReader(f):
                pid = pgs_id_of(r.get("accession", ""))
                if r.get("match_status") == "matched":
                    matched[pid] = matched.get(pid, 0) + int(float(r.get("count") or 0))
        score_dir = os.path.join(results, sampleset, "score")
        adjusted = os.path.join(score_dir, f"{sampleset}_pgs.txt.gz")
        plain = os.path.join(score_dir, "aggregated_scores.txt.gz")
        src = adjusted if os.path.isfile(adjusted) else plain
        if not os.path.isfile(src):
            raise Unreadable(f"pgsc_calc wrote no scores ({plain})")
        for r in read_gz_tsv(src):
            if r.get("sampleset") != sampleset:
                continue   # the reference panel's own samples
            pid = pgs_id_of(r.get("PGS", ""))
            sums[pid] = fmt_num(r.get("SUM"))
            if r.get("percentile_MostSimilarPop") not in (None, ""):
                v = fmt_num(r["percentile_MostSimilarPop"])
                pct[pid] = "NA" if v == "NA" else f"{float(v):.1f}"
        popsim = os.path.join(score_dir, f"{sampleset}_popsimilarity.txt.gz")
        if os.path.isfile(popsim):
            target = [r for r in read_gz_tsv(popsim)
                      if r.get("sampleset") == sampleset and r.get("REFERENCE", "False") != "True"]
            if len(target) != 1:
                raise Unreadable(f"{os.path.basename(popsim)} has {len(target)} rows for {sampleset}, expected 1")
            pops = target[0]
            for pid in sums:
                group[pid] = pops.get("MostSimilarPop") or "NA"
        if panel and chrx_scored:
            # chrX rows given to pgsc_calc that its matches did not keep:
            # "3 left out (not matched in the panel run)" for all of them,
            # "1 scored, 2 left out (not matched in the panel run)" for some
            xm = chrx_matched(results, sampleset)
            if xm is not None:
                for pid in chrx_scored:
                    n_x, got = int(chrx[pid].split()[0]), xm.get(pid, 0)
                    if not got:
                        chrx[pid] = f"{n_x} left out (not matched in the panel run)"
                    elif got < n_x:
                        chrx[pid] = f"{got} scored, {n_x - got} left out (not matched in the panel run)"
    rows = []
    for pid in sorted(totals):
        m, t = matched.get(pid, 0), totals[pid]
        a = applies.get(pid, "")
        check = "" if not a else (SEX_NOT_DECLARED if not sex else (SEX_APPLIES if sex == a else SEX_NOT_APPLICABLE))
        ev = evals.get(pid, [])
        # The effect is printed only for the group pgsc_calc placed the sample
        # in, and only when that group is one the score was evaluated in: an
        # odds ratio measured in one ancestry says little about another.
        g = group.get(pid, "NA")
        eff = next((f"{e} ({anc}, {ppm})" for anc, e, ppm in ev if anc == g), "")
        if check == SEX_NOT_APPLICABLE or pct.get(pid, "NA") == "NA" or below_log:
            eff = ""
        extra = [chrx.get(pid, ""), ambig.get(pid, ""), a, check, ";".join(x[0] for x in ev), eff, input_kind]
        if below_log:
            rows.append([labels.get(pid, pid), pid, "NA", "NA", str(t), rate.get(pid, "NA"), "NA", "NA"] + extra)
            continue
        # A score of the other sex keeps its raw sum and gets no percentile:
        # a percentile for an organ the person does not have reads as a risk.
        p = "NA" if check == SEX_NOT_APPLICABLE else pct.get(pid, "NA")
        rows.append([labels.get(pid, pid), pid, sums.get(pid, "NA"), str(m), str(t),
                     f"{100 * m / t:.1f}" if t else "0.0", p, group.get(pid, "NA")] + extra)
    ancestry = []
    if pops:
        pcs = sorted((k for k in pops if k.startswith("PC") and k[2:].isdigit()), key=lambda k: int(k[2:]))
        ancestry.append(("population", pops.get("MostSimilarPop", "")))
        ancestry.append(("population_low_confidence", pops.get("MostSimilarPop_LowConfidence", "")))
        for k in sorted(k for k in pops if k.startswith("RF_P_")):
            ancestry.append((f"probability_{k[5:]}", fmt_num(pops[k])))
        for k in pcs:
            ancestry.append((k, fmt_num(pops[k])))
    return rows, ancestry


def prs_table_main(argv):
    ap = argparse.ArgumentParser(prog="collect_summary.py prs-table",
                                 description="Write step 25's summary (and step 26's ancestry) table from pgsc_calc's output.")
    ap.add_argument("--sample", required=True)
    ap.add_argument("--results", required=True, help="pgsc_calc's --outdir")
    ap.add_argument("--sampleset", required=True, help="the sampleset name given to pgsc_calc")
    ap.add_argument("--scores", required=True, help="the folder prs-format wrote")
    ap.add_argument("--input-kind", required=True, choices=["gvcf", "vcf"],
                    help="gvcf: the score positions were genotyped from the gVCF; vcf: variant sites only")
    ap.add_argument("--panel", default="", help="name of the ancestry reference panel, when one was used")
    ap.add_argument("--zero-matches", action="store_true",
                    help="pgsc_calc stopped because no score variant is in the sample's genotypes")
    ap.add_argument("--below-threshold", metavar="LOG",
                    help="pgsc_calc stopped because every score matched under its minimum overlap; LOG is its console output")
    ap.add_argument("--sex", default="", choices=["", "male", "female"],
                    help="the declared sex: a score that applies to the other one gets no percentile")
    ap.add_argument("--out", required=True)
    ap.add_argument("--ancestry-out", help="written when pgsc_calc ran with the panel: PCs and population")
    a = ap.parse_args(argv)
    try:
        rows, ancestry = prs_table(a.results, a.sampleset, a.scores, a.input_kind, a.zero_matches, a.below_threshold,
                                   a.sex, bool(a.panel))
    except (Unreadable, OSError, ValueError, KeyError, EOFError, csv.Error, zlib.error) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        return 1
    # A score list without evaluations (an older or a custom one) gives the
    # table its older columns; the shipped list names some, so its table has
    # Evaluated_In and Effect before Input.
    evaluated = any(r[12] for r in rows)
    with open(a.out + ".tmp", "w", encoding="utf-8") as f:
        f.write("\t".join(PRS_COLUMNS[:-1] + PRS_EFFECT_COLUMNS + PRS_COLUMNS[-1:] if evaluated else PRS_COLUMNS) + "\n")
        for r in rows:
            f.write("\t".join(r if evaluated else r[:12] + r[14:]) + "\n")
    os.replace(a.out + ".tmp", a.out)
    if a.ancestry_out and ancestry:
        with open(a.ancestry_out + ".tmp", "w") as f:
            f.write("key\tvalue\n")
            for k, v in [("sample", a.sample), ("reference_panel", a.panel)] + ancestry:
                f.write(f"{k}\t{v}\n")
        os.replace(a.ancestry_out + ".tmp", a.ancestry_out)
    for r in rows:
        extra = f", percentile {r[6]} among {r[7]}" if r[6] != "NA" else ""
        if r[11] == SEX_NOT_APPLICABLE:
            extra = f", {r[10]}-specific: {SEX_NOT_APPLICABLE}, no percentile"
        elif r[11] == SEX_NOT_DECLARED:
            extra += f", {r[10]}-specific score: {SEX_NOT_DECLARED}"
        if r[13]:
            extra += f", effect in that group {r[13]}"
        if r[3] == "NA":
            print(f"  {r[0]} ({r[1]}): no sum, under pgsc_calc's minimum overlap ({r[5]}% of {r[4]} variants matched)")
        else:
            print(f"  {r[0]} ({r[1]}): sum {r[2]}, {r[3]} of {r[4]} variants matched ({r[5]}%){extra}")
    return 0


def count_vcf(p):
    total = passed = 0
    for r in vcf_records(p):
        total += 1
        if len(r) > 6 and r[6] == "PASS":
            passed += 1
    return total, passed


def sec_ancestry_fine(d, s):
    """Step 26b: <s>_ancestry_fine.tsv (level, code, group, fraction, se, low,
    high; levels continent, region and painting, plus left_out for a group the
    model leaves out and not_painted for a chromosome RFMix failed on) and
    the painting picture."""
    p = first_existing(d, [f"ancestry_fine/{s}_ancestry_fine.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    need = {"level", "code", "group", "fraction", "se", "low", "high"}
    if not rows or not need <= set(rows[0]):
        raise Unreadable(f"{os.path.basename(p)} has no rows or lacks a column of {sorted(need)}")
    out = {"continent": [], "region": [], "painting": [], "left_out": [], "not_painted": [],
           "min_references": ancestry_fine_report.MIN_REFERENCES}
    for r in rows:
        if r["level"] in out:
            out[r["level"]].append({"code": r["code"], "group": r["group"], "fraction": num(r["fraction"]),
                                    "low": num(r["low"]), "high": num(r["high"])})
    # The painting picture travels in the summary, so the HTML report shows it
    # from the summary alone (render_report.py --summary).
    png = first_existing(d, [f"ancestry_fine/{s}_painting.png"])
    out["painting_png"] = None
    out["legend"] = [{"code": c, "group": ancestry_fine_report.label("continent", c),
                      "colour": "#%02x%02x%02x" % rgb} for c, rgb in ancestry_fine_report.COLOURS.items()]
    if png and os.path.getsize(png) < 2_000_000:
        with open(png, "rb") as f:
            out["painting_png"] = base64.b64encode(f.read()).decode()
    return p, out


def sec_ancient(d, s):
    """Step 44: ancient/<s>_ancient.tsv (the nearest ancient groups on the AADR
    panel's PCs), <s>_ancient_info.tsv (release, sites used, people, the
    sample's PC1 and PC2) and the plot, which travels in the summary so the
    HTML report shows it from the summary alone."""
    p = first_existing(d, [f"ancient/{s}_ancient.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    need = {"rank", "group", "n", "distance", "dates", "country", "locality"}
    if not rows or not need <= set(rows[0]):
        raise Unreadable(f"{os.path.basename(p)} has no rows or lacks a column of {sorted(need)}")
    info = {}
    ip = first_existing(d, [f"ancient/{s}_ancient_info.tsv"])
    if ip:
        info = {r.get("key"): r.get("value") for r in read_tsv(ip)}
    out = {"nearest": [{k: r[k] for k in ("rank", "group", "n", "distance", "dates", "country", "locality")} for r in rows],
           "release": info.get("release") or ".", "sites_used": info.get("sites_used") or ".",
           "panel_sites": info.get("panel_sites") or ".", "ancient_people": info.get("ancient_people") or ".",
           "ancient_groups": info.get("ancient_groups") or ".", "present_day_people": info.get("present_day_people") or ".",
           "plot_svg": None}
    svg = first_existing(d, [f"ancient/{s}_ancient_pca.svg"])
    if svg and os.path.getsize(svg) < 3_000_000:
        with open(svg, "rb") as f:
            out["plot_svg"] = base64.b64encode(f.read()).decode()
    return p, out


def sec_manta(d, s):
    p = first_existing(d, ["manta/results/variants/diploidSV.vcf.gz", "manta/*diploidSV.vcf.gz",
                           "manta/**/diploidSV.vcf.gz"])
    if not p:
        return None, {}
    t, ps = count_vcf(p)
    return p, {"total": t, "pass": ps}


def sec_delly(d, s):
    p = first_existing(d, [f"delly/{s}_sv.vcf.gz"])
    if not p:
        return None, {}
    t, ps = count_vcf(p)
    return p, {"total": t, "pass": ps}


# A CNVpytor call is kept in the filtered count when fewer than half of its
# reads have mapping quality 0 (q0, column 9) and fewer than half of its bases
# are N (pN, column 10), it touches no centromere or heterochromatin band
# (MASKED_BANDS), and it is not on chrY of a declared female (collect() takes
# those out, once the declared sex is known).
CNV_MAX_Q0 = CNV_MAX_PN = 0.5


def sec_cnvpytor(d, s):
    p = first_existing(d, [f"cnvpytor/{s}_cnvs.txt"])
    if not p:
        return None, {}
    bands = masked_bands()
    total = dels = dups = sig = kept = kept_y = 0
    dropped = {"q0": 0, "gap": 0, "band": 0}
    with open(p, errors="replace") as f:
        for line in f:
            c = line.split()
            if not c:
                continue
            total += 1
            dels += c[0] == "deletion"
            dups += c[0] == "duplication"
            m = re.fullmatch(r"(\S+):(\d+)-(\d+)", c[1]) if len(c) > 4 else None
            try:
                if float(c[4]) >= 0.01:
                    continue
            except (IndexError, ValueError):
                continue
            sig += 1
            q0, pn = (num(c[8]), num(c[9])) if len(c) > 9 else (None, None)
            if q0 is not None and q0 >= CNV_MAX_Q0:
                dropped["q0"] += 1
            elif pn is not None and pn >= CNV_MAX_PN:
                dropped["gap"] += 1
            elif m and masked_overlap(bands, m.group(1), int(m.group(2)), int(m.group(3))):
                dropped["band"] += 1
            else:
                kept += 1
                kept_y += bool(m) and m.group(1) in ("chrY", "Y")
    return p, {"total": total, "deletions": dels, "duplications": dups, "significant": sig,
               "significant_filtered": kept, "filtered_chry": kept_y, "dropped": dropped}


def cnvpytor_sex_filter(data, declared_sex):
    """Take the chrY calls out of the filtered count of a declared female;
    with no declared sex they stay, and the data says so."""
    if "significant_filtered" not in data:
        return
    data["declared_sex"] = declared_sex or None
    if declared_sex == "female":
        data["dropped"]["chry_female"] = data["filtered_chry"]
        data["significant_filtered"] -= data["filtered_chry"]


def sec_sv_consensus(d, s):
    p = first_existing(d, [f"sv_merged/{s}_sv_consensus.vcf.gz"])
    if not p:
        return None, {}
    t, _ = count_vcf(p)
    return p, {"consensus": t}


def sec_annotsv(d, s):
    """Step 05's class 4 and 5 SVs (annotsv/<s>_sv_class4_5_checked.tsv,
    bin/annotsv_postfilter.py), each with how many callers support it from
    step 22's consensus (sv_merged/<s>_sv_consensus.vcf.gz), joined here
    because step 22 runs after step 05 (and beside ANNOTSV in Nextflow)."""
    p = first_existing(d, [f"annotsv/{s}_sv_class4_5_checked.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p) if os.path.getsize(p) else []
    need = {"SV_chrom", "SV_start", "SV_end", "SV_type", "ACMG_class", "Caller", "Exon_edge", "Population_AF", "Common"}
    if rows and not need <= set(rows[0]):
        raise Unreadable(f"{os.path.basename(p)} lacks a column of {sorted(need)}")
    cons = first_existing(d, [f"sv_merged/{s}_sv_consensus.vcf.gz"])
    cons_note = None
    try:
        support = annotsv_postfilter.consensus_support(rows, cons)
    except (OSError, EOFError, zlib.error, UnicodeDecodeError) as e:
        support = annotsv_postfilter.consensus_support(rows, None)
        cons_note = f"{os.path.basename(cons)} could not be read: {e}"
        cons = None
    out = []
    for r, sp in zip(rows, support):
        try:
            af = float(r["Population_AF"]) if r["Population_AF"] else None
        except ValueError:
            af = None
        genes = [g for g in (r.get("Gene_name") or "").split(";") if g]
        fmt = annotsv_postfilter.sample_format(r)
        out.append({"sv": f"{r['SV_chrom']}:{r['SV_start']}-{r['SV_end']}", "type": r["SV_type"],
                    # the call's genotype and phase set (PS), for the two-hit rule (two_hits)
                    "gt": fmt.get("GT") or None, "ps": fmt.get("PS") or None,
                    "length": r.get("SV_length") or "", "acmg_class": r["ACMG_class"], "genes": genes,
                    "caller": r["Caller"], "support": sp["Support"],
                    "support_callers": sp["Support_callers"].split(",") if sp["Support_callers"] else [],
                    "exon_edge": r["Exon_edge"], "exon_edge_detail": r.get("Exon_edge_detail") or "",
                    "population_af": af, "population_source": r.get("Population_source") or "",
                    "common": r["Common"] == "yes", "annotsv_class": r.get("AnnotSV_class") or r["ACMG_class"],
                    "rule": r.get("Rule") or "", "rule_detail": r.get("Rule_detail") or ""})
    # A call step 05's rules moved from AnnotSV's class 4 or 5 to class 3
    # (bin/annotsv_postfilter.py: a coding tandem repeat, an IG/TR locus, the
    # MHC, a pseudogene alias, no read-depth support) is listed apart with
    # the reason, and the class 4-5 counts leave it out.
    moved = [x for x in out if str(x["acmg_class"]) not in ("4", "5")]
    out = [x for x in out if str(x["acmg_class"]) in ("4", "5")]
    return p, {"count": len(out), "rows": out, "downgraded": moved,
               "consensus": os.path.relpath(cons, d) if cons else None, "consensus_note": cons_note,
               "single_caller": sum(x["support"] == 1 for x in out),
               "exon_edge": sum(x["exon_edge"] == "yes" for x in out),
               "exon_edge_checked": not any(x["exon_edge"] == "not checked" for x in out),
               "common": sum(x["common"] for x in out)}


STR_RANGES_TSV = os.path.join(os.path.dirname(os.path.abspath(__file__)), "str_intermediate_ranges.tsv")
# How well ExpansionHunter sized an allele, from its FORMAT SO, ADSP and
# REPCI fields, worst first; a locus takes the worst of its alleles.
STR_SIZING = ["no_call", "flanking_only", "in_repeat", "low_confidence", "sized", "unknown"]
STR_MIN_SPANNING = 2


def str_ranges(path=STR_RANGES_TSV):
    """bin/str_intermediate_ranges.tsv by locus (tests/test_str_ranges.py checks it)."""
    rows, head = {}, None
    with open(path) as f:
        for line in f:
            if line.startswith("#") or not line.strip():
                continue
            c = line.rstrip("\n").split("\t")
            if head is None:
                head = c
                continue
            r = dict(zip(head, c))
            rows[r["locus"]] = r
    return rows


def str_int(x):
    try:
        return int(x)
    except (TypeError, ValueError):
        return None


def str_allele_sizing(so, adsp, ci, cn, path_min):
    """One allele's sizing: sized, low_confidence (fewer than
    STR_MIN_SPANNING spanning reads, or a confidence interval that reaches
    the catalog's pathogenic minimum from below it), in_repeat (in-repeat
    reads only: longer than a read), flanking_only, no_call, or unknown
    (no SO field)."""
    if str_int(cn) is None:
        return "no_call"
    if not so or so == ".":
        return "unknown"
    if so == "FLANKING":
        return "flanking_only"
    if so == "INREPEAT":
        return "in_repeat"
    hi = str_int(ci.split("-")[-1]) if ci else None
    if (str_int(adsp) is not None and str_int(adsp) < STR_MIN_SPANNING) or (
            path_min is not None and hi is not None and str_int(cn) < path_min <= hi):
        return "low_confidence"
    return "sized"


def str_field(fmt, val, key):
    return val[fmt.index(key)] if key in fmt and fmt.index(key) < len(val) else ""


def str_split(x, n):
    """A per-allele FORMAT value as n parts; one value serves every allele."""
    parts = x.split("/") if x else []
    return parts if len(parts) == n else (parts * n if len(parts) == 1 else [""] * n)


def rfc1_motif_row(d, s, after):
    """Step 9b's RFC1 motif count (bin/rfc1_motif.py), unless it is older than
    the Stranger file it was made from."""
    p = first_existing(d, [f"expansion_hunter/{s}_rfc1_motif.tsv"])
    if not p or (after and os.path.getmtime(p) < os.path.getmtime(after)):
        return None
    rows = read_tsv(p)
    if not rows:
        return None
    r = rows[0]
    return {k: r.get(k, ".") for k in ("verdict", "note", "dominant_motif", "dominant_class",
                                       "dominant_share", "motif_shares", "repeat_reads")}


def str_motif_rows(d, s, after):
    """Step 9b's motif count of the other flagged motif loci (the FAME loci,
    SCA31, SCA37; expansion_hunter/<id>_str_motifs.tsv, bin/rfc1_motif.py
    --flagged) by locus, unless the table is older than the Stranger file it
    was made from."""
    p = first_existing(d, [f"expansion_hunter/{s}_str_motifs.tsv"])
    if not p or (after and os.path.getmtime(p) < os.path.getmtime(after)):
        return {}
    return {r["locus"]: {k: r.get(k, ".") for k in ("verdict", "note", "dominant_motif", "dominant_class",
                                                    "dominant_share", "motif_shares", "repeat_reads")}
            for r in read_tsv(p) if r.get("locus")}


def str_not_graded(d, s):
    """{locus: reason} of the loci the catalog ExpansionHunter ran on gives no
    thresholds (its NoThresholds field, bin/strchive_convert.py: MUC1, CEL,
    TYMS and PLIN4 among them), from the copy step 09 keeps beside its VCF
    (expansion_hunter/<id>_eh_catalog.json). Empty without the copy, or with
    one that is not a catalog (Stranger's bundled catalog has no such loci)."""
    p = first_existing(d, [f"expansion_hunter/{s}_eh_catalog.json"])
    if not p:
        return {}
    try:
        with open(p, encoding="utf-8") as f:
            cat = json.load(f)
    except (OSError, ValueError):
        return {}
    if not isinstance(cat, list):
        return {}
    return {e["LocusId"]: str(e["NoThresholds"]) for e in cat
            if isinstance(e, dict) and e.get("LocusId") and e.get("NoThresholds")}


# The repeats STRchive leaves out (assets/str_extra_loci.json), read from
# their counts when step 09 sized them (STR_EXTRA_LOCI=yes; docs/09-str-expansions.md).
STR_EXTRA_SETUP = "STR_EXTRA_LOCI=yes ./scripts/setup.sh --strchive-data <genome_dir>"
# What a section state other than ok says about a table another step wrote.
OTHER_STEP_STATE = {"stale": "from an earlier run", "failed": "its step failed in the last run",
                    "unreadable": "unreadable"}
# CFTR intron 9 (legacy IVS8): ExpansionHunter writes the TG and the poly-T
# tract as two records, so which TG count sits on the copy of which T count
# is known only when one of the two is the same on both copies. R117H
# (rs78655421, chr7:117530975 G>A) is step 42's traits row, 17.6 kb from the
# tract, farther than a read pair spans: with both heterozygous its phase is
# not known. Deignan et al. 2020 (PMID 32404922).
CFTR_R117H = "rs78655421"
CFTR_SOURCE = "Deignan et al. 2020, PMID 32404922"
# UGT1A1 promoter: ExpansionHunter counts the TA units from the first T of
# A(TA)6TAA, so its count is one more than the allele's TA number.
UGT1A1_ALLELES = {6: "*36", 7: "*1", 8: "*28", 9: "*37"}
UGT1A1_DECREASED = ("*28", "*37")
UGT1A1_SOURCE = "Gammal et al. 2016, PMID 26417955"
HTTLPR_READING = ("not tested by any step. CPIC's 2023 SSRI guideline (Bousman et al. 2023, PMID 37032427) found "
                  "that the data on SLC6A4 do not support its use in antidepressant prescribing, so a result would "
                  "carry no recommendation.")


def str_counts(item):
    """The two repeat counts of a sized record, sorted, or None."""
    if not item or item["sizing"] not in ("sized", "unknown"):
        return None
    c = [str_int(x) for x in item["repeat_count"].split("/")]
    if len(c) == 1 and c[0] is not None:
        c = c * 2  # a one-copy call (a haploid contig): not expected on chr7 or chr2
    return sorted(c) if len(c) == 2 and None not in c else None


def trait_genotype(d, s, rsid):
    """(status, genotype) of a traits table row of step 42 (traits/<id>_traits.tsv), or None."""
    p = first_existing(d, [f"traits/{s}_traits.tsv"])
    if not p:
        return None
    try:
        row = next((r for r in read_tsv(p) if r.get("rsid") == rsid), None)
    except (OSError, csv.Error):
        return None
    return (row.get("status") or "", row.get("genotype") or "") if row else None


def r117h_genotype(d, s, state):
    """(genotype or None, why not read, note) of R117H in step 42's traits
    table, given the traits section's state: a table of a failed or
    unreadable step is not read, one from an earlier run is read and said so."""
    if state in ("failed", "unreadable"):
        return None, f"step 42's traits table is {OTHER_STEP_STATE[state]}", ""
    r = trait_genotype(d, s, CFTR_R117H)
    if not r or r[0] != "called":
        return None, "no step 42 traits row for rs78655421 was called", ""
    return r[1], "", (" (R117H from step 42's traits table of an earlier run)" if state == "stale" else "")


def cftr_reading(by, d, s, traits_state="ok"):
    tg, t = by.get("CFTR_TG"), by.get("CFTR_T")
    out = {"locus": "CFTR", "repeat": "intron 9 (legacy IVS8) TG and poly-T tract", "source": CFTR_SOURCE}
    if not tg or not t:
        return dict(out, status="not_sized", genotype=".", phase=".",
                    reading=f"not in step 09's output: {STR_EXTRA_SETUP} adds the tract to the catalog "
                            "(docs/09-str-expansions.md).")
    a, b = str_counts(tg), str_counts(t)
    if not a or not b:
        return dict(out, status="not_sized", genotype=".", phase=".",
                    reading=f"ExpansionHunter did not size both tracts (TG {tg['repeat_count']}, {tg['sizing']}; "
                            f"poly-T {t['repeat_count']}, {t['sizing']}): not read.")
    if a[0] == a[1] or b[0] == b[1]:
        pairs = [f"TG{a[0]}-{b[0]}T/TG{a[1]}-{b[1]}T"]
    else:
        pairs = [f"TG{a[0]}-{b[0]}T/TG{a[1]}-{b[1]}T", f"TG{a[1]}-{b[0]}T/TG{a[0]}-{b[1]}T"]
    genotype = pairs[0] if len(pairs) == 1 else f"TG{a[0]}/TG{a[1]} {b[0]}T/{b[1]}T"
    out.update(status="read", genotype=genotype, pairings=pairs, tg=a, poly_t=b)
    five = b.count(5)
    other = sorted({n for n in b if n not in (5, 7, 9)})
    odd = (f" A poly-T length of {', '.join(f'{n}T' for n in other)} is none of the usual 5T, 7T and 9T."
           if other else "")
    r117, why, note = r117h_genotype(d, s, traits_state)
    if not five:
        # No 5T risk of its own, but the poly-T still sets how R117H reads:
        # with no 5T copy, an R117H copy sits on a 7T or 9T (or other) copy.
        line = (f"No 5T copy ({b[0]}T/{b[1]}T): the tract adds no 5T allele of its own. The poly-T length still "
                "matters beside R117H, which is milder and of variable penetrance on a 7T or 9T copy and CF-causing "
                "with 5T in cis.")
        if r117 is None:
            line += f" R117H was not read ({why})."
        elif r117.split("/").count("A") == 2:
            copies = f"both {b[0]}T copies" if b[0] == b[1] else f"both copies ({b[0]}T and {b[1]}T)"
            line += f" R117H ({r117}) is read: it sits on {copies}.{note}"
        elif r117.split("/").count("A"):
            copies = f"{b[0]}T" if b[0] == b[1] else f"{b[0]}T or {b[1]}T"
            line += f" R117H ({r117}) is read: with no 5T copy it sits on a {copies} copy.{note}"
        else:
            line += f" No R117H copy.{note}"
        return dict(out, five_t=False, phase="not_needed", r117h=r117 or ".", reading=line + odd)
    if len(pairs) == 1:
        # One TG count on both copies, or 5T on both copies with their two TG counts.
        tgs = [a[0]] if a[0] == a[1] else a
        tg_text = " and ".join(f"TG{x}" for x in tgs) + (" on the 5T copy" if five == 1 else " on the 5T copies")
    else:
        tg_text = f"TG{a[0]} or TG{a[1]} on the 5T copy, not told by the counts ({' or '.join(pairs)})"
    line = (f"5T: {five} copy{'s' if five > 1 else ''} with the 5T poly-T tract, {tg_text}. 5T lowers the splicing "
            "of exon 10 and is a variable-penetrance allele: with TG12 or TG13 on the same copy it more often "
            "causes CBAVD or another CFTR-related disorder when the other copy carries a CF-causing variant.")
    if r117 is None:
        phase = "not_read"
        text = f" R117H was not read ({why}), so whether it sits on the 5T copy is not known."
    else:
        n = r117.split("/").count("A")
        if n == 0:
            phase = "not_needed"
            text = " No R117H copy: the 5T is read on its own."
        elif n == 2 or five == 2:
            phase = "in_cis"
            text = (f" R117H ({r117}) and 5T share a copy: R117H with 5T in cis is a CF-causing allele of "
                    "variable penetrance.")
        else:
            phase = "unknown"
            text = (" One R117H copy and one 5T copy: phase unknown. R117H lies 17.6 kb from the tract, farther "
                    "than a read pair spans, so whether they share a copy (in cis, a CF-causing allele of variable "
                    "penetrance) or not (in trans) is not read; testing a parent tells them apart.")
        text += note
    return dict(out, five_t=True, phase=phase, r117h=r117 or ".", reading=line + text + odd)


def ugt1a1_reading(by, d, s, cpic_state="ok"):
    item = by.get("UGT1A1")
    out = {"locus": "UGT1A1", "repeat": "promoter TA repeat", "source": UGT1A1_SOURCE}
    if not item:
        return dict(out, status="not_sized", genotype=".", pharmcat=".", agreement="not_compared",
                    reading=f"not in step 09's output: {STR_EXTRA_SETUP} adds the repeat to the catalog "
                            "(docs/09-str-expansions.md). PharmCAT (step 7) still reads UGT1A1 *28.")
    c = str_counts(item)
    if not c:
        return dict(out, status="not_sized", genotype=".", pharmcat=".", agreement="not_compared",
                    reading=f"ExpansionHunter did not size the repeat ({item['repeat_count']}, {item['sizing']}): "
                            "not read.")
    names = [UGT1A1_ALLELES.get(n, f"TA{n - 1}") for n in c]
    genotype = "/".join(names)
    dec = sum(n in UGT1A1_DECREASED for n in names)
    ta = "/".join(f"TA{n - 1}" for n in c)
    out.update(status="read", genotype=genotype, ta_units=c)
    reading = f"{genotype} ({ta}, {c[0]}/{c[1]} TA units as ExpansionHunter counts them): "
    if not all(n in UGT1A1_ALLELES for n in c):
        # A length with no named allele has no function to read, and so no
        # count of decreased-function copies to set beside PharmCAT's.
        return dict(out, pharmcat=".", agreement="not_compared",
                    reading=reading + "a TA length outside TA5 to TA8 is not a named allele here, so its function "
                                      "is not read and it is not compared with PharmCAT (step 7).")
    reading += ("no decreased-function TA allele." if not dec else
                f"{dec} decreased-function allele{'s' if dec > 1 else ''} (*28 and *37 are decreased function).")
    # PharmCAT's UGT1A1 call (step 27's table, with the reference-block fill
    # of bin/pgx_refblock_fill.py): its count of *28 and *37 copies.
    pc, agreement = ".", "not_compared"
    p = first_existing(d, [f"cpic/{s}_phenotypes.tsv"]) if cpic_state not in ("failed", "unreadable") else None
    row = None
    if p:
        try:
            row = next((r for r in read_tsv(p) if r.get("Gene") == "UGT1A1"), None)
        except (OSError, csv.Error):
            row = None
    if row:
        pc = (row.get("Diplotype") or "").strip() or "."
        alleles = pc.split("/")
        if len(alleles) == 2 and "possible" not in pc and (row.get("Status") or "") != "not called":
            # A combination allele reads "*80+*28" or "[*80 + *28]".
            n = sum(any(x in re.split(r"[+ \[\]]", a) for x in UGT1A1_DECREASED) for a in alleles)
            agreement = "agrees" if n == dec else "differs"
    old = " (step 27's table is from an earlier run)" if cpic_state == "stale" and agreement != "not_compared" else ""
    if agreement == "agrees":
        reading += f" PharmCAT (step 7) called {pc}: the same number of *28 and *37 copies{old}."
    elif agreement == "differs":
        reading += (f" PharmCAT (step 7) called {pc}, a different number of *28 and *37 copies{old}: read the reads "
                    "at chr2:233760234 before using either.")
    elif cpic_state in ("failed", "unreadable"):
        reading += f" Not compared with PharmCAT (step 27's table is {OTHER_STEP_STATE[cpic_state]})."
    else:
        reading += " Not compared with PharmCAT (no unambiguous UGT1A1 call from step 27)."
    return dict(out, pharmcat=pc, agreement=agreement, reading=reading)


def str_readings(by, d, s, states=None):
    """What the counts of the extra repeats mean, one entry per locus, and the
    5-HTTLPR line that no step can test. states: the traits and cpic section
    states, once collect() knows them (a stale, failed or unreadable table is
    said so, not read as current)."""
    states = states or {}
    return [cftr_reading(by, d, s, states.get("traits", "ok")), ugt1a1_reading(by, d, s, states.get("cpic", "ok")),
            {"locus": "SLC6A4", "repeat": "5-HTTLPR promoter length polymorphism", "status": "not_tested",
             "genotype": ".", "reading": HTTLPR_READING, "source": "Bousman et al. 2023, PMID 37032427"}]


def sec_expansions(d, s):
    eh = first_existing(d, [f"expansion_hunter/{s}_eh.vcf", f"expansion_hunter/{s}.vcf",
                            "expansion_hunter/*_eh.vcf"])
    st = first_existing(d, [f"expansion_hunter/{s}_eh_stranger.vcf", "expansion_hunter/*_eh_stranger.vcf"])
    # Stranger (step 9b) copies ExpansionHunter's records and adds STR_STATUS.
    # Its file is read unless it is older than the ExpansionHunter file it
    # would have been made from; then the report says to rerun step 9b.
    outdated = bool(st and eh and os.path.getmtime(st) < os.path.getmtime(eh))
    if outdated:
        st = None
    p = st or eh
    if not p:
        return None, {}
    ranges = str_ranges()
    rfc1 = rfc1_motif_row(d, s, st) if st else None
    motifs = str_motif_rows(d, s, st) if st else {}
    ungraded = str_not_graded(d, s) if st else {}
    loci, tested, flagged, no_status, not_graded = [], 0, [], 0, []
    for r in vcf_records(p):
        if len(r) < 10:
            continue
        tested += 1
        info = info_map(r[7])
        # A locus of several repeats (CFTR: CFTR_TG and CFTR_T) writes one
        # record per repeat, each with the locus as REPID and its own VARID:
        # the record is named by its VARID and looked up by that name in the
        # ranges and motif tables (the pathogenic repeat of a multi-repeat
        # locus carries the locus name), and in the catalog's NoThresholds by
        # its locus (REPID). A one-repeat locus has VARID = REPID.
        lid = info.get("REPID") or info.get("VARID") or f"{r[0]}:{r[1]}"
        rep = info.get("VARID") or lid
        fmt = r[8].split(":")
        val = r[9].split(":")
        repcn = str_field(fmt, val, "REPCN")
        counts = repcn.split("/") if repcn and repcn != "." else []
        n = max(len(counts), 1)
        path_min = str_int(info.get("STR_PATHOLOGIC_MIN"))
        sizing = [str_allele_sizing(so, ad, ci, cn, path_min) for so, ad, ci, cn in
                  zip(str_split(str_field(fmt, val, "SO"), n), str_split(str_field(fmt, val, "ADSP"), n),
                      str_split(str_field(fmt, val, "REPCI"), n), counts or ["."])]
        status = [x for x in info.get("STR_STATUS", "").split(",") if x and x != "."] if st else []
        # Stranger writes '.' where its catalog has no value.
        item = {"locus": rep, "repeat_count": repcn or ".", "status": ",".join(status),
                "normal_max": info.get("STR_NORMAL_MAX", "").strip("."),
                "pathologic_min": info.get("STR_PATHOLOGIC_MIN", "").strip("."),
                "support": str_field(fmt, val, "SO"), "adsp": str_field(fmt, val, "ADSP"),
                "ci": str_field(fmt, val, "REPCI"), "sizing": min(sizing, key=STR_SIZING.index),
                "allele_sizing": "/".join(sizing), "zone": "", "zone_label": "", "zone_citation": ""}
        rng = ranges.get(rep)
        sized = [str_int(c) for c, z in zip(counts, sizing) if z in ("sized", "unknown") and str_int(c) is not None]
        if rng and any(int(rng["intermediate_min"]) <= c <= int(rng["intermediate_max"]) for c in sized):
            item.update(zone=f"{rng['intermediate_min']}-{rng['intermediate_max']}", zone_label=rng["label"],
                        zone_citation=rng["citation"])
        if rep == "RFC1" and st and any(x != "normal" for x in status):
            item["rfc1_motif"] = rfc1 or {"verdict": "not_checked",
                                          "note": "no motif count beside Stranger's output (step 9b runs it "
                                                  "with the sample's BAM): confirm the flag with repeat-primed PCR"}
        if rep != "RFC1" and rep in motifs and any(x != "normal" for x in status):
            item["motif"] = motifs[rep]
        if st and not status and lid in ungraded:
            item["not_graded"] = ungraded[lid]
        loci.append(item)
        if st:
            # One value per record, the most severe of normal, pre_mutation and
            # full_mutation; a locus outside Stranger's catalog has none.
            if not status and "not_graded" in item:
                if lid not in not_graded:
                    not_graded.append(lid)
            elif not status:
                no_status += 1
            elif any(x != "normal" for x in status):
                f = {"locus": item["locus"], "repeat_count": item["repeat_count"], "status": item["status"],
                     "sizing": item["sizing"]}
                if "rfc1_motif" in item:
                    f["rfc1_motif"] = item["rfc1_motif"]
                if "motif" in item:
                    f["motif"] = item["motif"]
                flagged.append(f)
    by = {l["locus"]: l for l in loci}
    # Every locus, the five key loci first, then by name.
    loci.sort(key=lambda l: (EH_LOCI.index(l["locus"]) if l["locus"] in EH_LOCI else len(EH_LOCI), l["locus"]))
    out = {"records": tested, "key_loci": [{"locus": k, "repeat_count": by[k]["repeat_count"] if k in by
                                            else "not in output"} for k in EH_LOCI],
           "loci": loci, "stranger": bool(st), "stranger_outdated": outdated,
           "unsized": [l["locus"] for l in loci if l["sizing"] not in ("sized", "unknown")],
           "intermediate": [l["locus"] for l in loci if l["zone"]],
           "readings": str_readings(by, d, s)}
    if st:
        # not_graded: {locus: why} of the loci the catalog gives no thresholds,
        # so Stranger gives them no status: not graded, never normal.
        out.update(flagged=flagged, no_status=no_status, not_graded={k: ungraded[k] for k in not_graded})
    return p, out


# TelomereHunter's estimate rests on its intratelomeric reads. Below this many
# per 30x of depth, about 30 times fewer than a typical 30x short-read genome
# has, the value is too thin to compare (a library or platform problem, or
# reads too short for its repeat threshold). The floor scales with the depth
# TelomereHunter itself counted (total_reads x read length / GENOME_BASES).
TELOMERE_FLOOR_PER_30X = 2000
GENOME_BASES = 3.1e9


def sec_telomere(d, s):
    p = first_existing(d, [f"telomere/{s}/{s}/{s}_summary.tsv", f"telomere/{s}/{s}_summary.tsv",
                           "telomere/**/*_summary.tsv"])
    if not p:
        return None, {}
    with open(p, errors="replace") as f:
        lines = [l.rstrip("\n").split("\t") for l in f if l.strip()]
    if len(lines) < 2:
        return p, {"tel_content": "."}
    head, row = lines[0], lines[1]
    # By header name; older TelomereHunter versions put it in column 11.
    i = head.index("tel_content") if "tel_content" in head else 10
    out = {"tel_content": row[i] if i < len(row) else "."}

    def col(*names):
        for n in names:
            if n in head and head.index(n) < len(row):
                return row[head.index(n)]
        return None
    # TelomereHunter 1.1.0 names it intratel_reads.
    intratel = num(col("intratel_reads", "intratelomeric_reads"))
    total = num(col("total_reads"))
    lengths = [num(x) for x in re.split(r"[,;| ]+", col("read_lengths", "read_length") or "")]
    lengths = [x for x in lengths if x]
    if intratel is None:
        return p, dict(out, read_check="no intratelomeric read count in the summary")
    out["intratelomeric_reads"] = int(intratel)
    if not total or not lengths:
        return p, dict(out, read_check="depth unknown (no total_reads or read_lengths in the summary): not checked")
    depth = total * max(lengths) / GENOME_BASES
    floor = TELOMERE_FLOOR_PER_30X * depth / 30
    out.update(depth=round(depth, 1), read_floor=int(round(floor)), low_reads=intratel < floor)
    return p, out


# GRCh38's centromeres and heterochromatin: the acen, gvar and stalk rows of
# UCSC's cytoBand.hg38.txt (the file setup.sh installs as
# reference/cytoBand.hg38.txt, sha256 e514e48b... in scripts/lib/common.sh),
# neighbouring rows merged, as chrom:start-end with a 0-based start and an
# exclusive end. Reads there map badly, so a run of homozygosity or a copy
# number call over them is mostly an artefact. Regenerate with
#   gzip -cd cytoBand.txt.gz | awk '$1 ~ /^chr([0-9]+|X|Y)$/ && $5 ~ /^(acen|gvar|stalk)$/' \
#     | sort -k1,1V -k2,2n | awk '{ if ($1 == c && $2 == e) e = $3; else { if (c != "") printf "%s:%d-%d ", c, s, e;
#       c = $1; s = $2; e = $3 } } END { print c ":" s "-" e }'
# scripts/11-roh-analysis.sh and modules/local/roh/main.nf carry the same
# line (tests/test_collect_summary.py checks that the three agree).
MASKED_BANDS = ("chr1:121700000-143200000 chr2:91800000-96000000 chr3:87800000-98600000 chr4:48200000-51800000 "
                "chr5:46100000-51400000 chr6:58500000-62600000 chr7:58100000-62100000 chr8:43200000-47200000 "
                "chr9:42200000-61500000 chr10:38000000-41600000 chr11:51000000-55800000 chr12:33200000-37800000 "
                "chr13:0-18900000 chr14:0-18200000 chr15:0-20500000 chr16:35300000-47000000 "
                "chr17:22700000-27400000 chr18:15400000-21500000 chr19:19900000-31900000 chr20:25700000-30400000 "
                "chr21:0-13000000 chr22:0-17400000 chrX:58100000-63800000 chrY:10300000-10600000 "
                "chrY:26600000-57227415")


def masked_bands(text=MASKED_BANDS):
    """{chrom: [(start, end), ...]} of MASKED_BANDS, 0-based, end exclusive."""
    out = {}
    for x in text.split():
        c, _, se = x.partition(":")
        a, _, b = se.partition("-")
        out.setdefault(c, []).append((int(a), int(b)))
    return out


def masked_overlap(bands, chrom, start, end):
    """Bases of the 1-based, inclusive span start-end inside the masked bands."""
    return sum(max(0, min(end, b) - max(start - 1, a)) for a, b in bands.get(chrom, ()))


# GRCh38's chr1 to chr22 (UCSC hg38.chrom.sizes), and the part of them outside
# MASKED_BANDS: the length F_ROH is a share of, since a run's part in a band is
# not counted either.
AUTOSOME_BP = 2875001522
ROH_MIN_BP = 5e6     # the segments F_ROH and the one-chromosome flag count (docs/11-roh-analysis.md)
# One chromosome carrying most of the long runs: at least UPD_MIN_BP of it in
# segments of 5 Mb or more, and more than UPD_SHARE of the length of all such
# autosomal runs.
# Our thresholds: a run over 20 Mb on one chromosome is the long-standing sign
# of isodisomy, and related parents spread their runs over many chromosomes.
UPD_MIN_BP = 20e6
UPD_SHARE = 0.5
# The rows of docs/11-roh-analysis.md's table: F of a child of the parents.
# A value is placed at the row whose F is nearest on a log scale (the
# boundaries are the geometric means of neighbouring rows), and under 1/128
# (half the second-cousin F) as no sign of related parents.
ROH_BANDS = [(1 / 128, "no sign of related parents"),
             ((1 / 64 * 1 / 32) ** 0.5, "second cousins (F 1/64)"),
             ((1 / 32 * 1 / 16) ** 0.5, "first cousins once removed (F 1/32)"),
             ((1 / 16 * 1 / 8) ** 0.5, "first cousins (F 1/16)"),
             ((1 / 8 * 1 / 4) ** 0.5, "half siblings, uncle and niece, double first cousins (F 1/8)"),
             (float("inf"), "parent and child, full siblings (F 1/4)")]

# assets/imprinted_regions.tsv, the same rows (tests/test_collect_summary.py
# checks that the two agree). Held here too because the Nextflow report task
# sees bin/ but not assets/. (chrom, start, end, region, maternal UPD,
# paternal UPD, PMIDs); start and end None for the whole chromosome.
IMPRINTED_REGIONS = [
    ("chr6", None, None, "6q24 (PLAGL1, HYMAI)", "-", "Transient neonatal diabetes mellitus (6q24-related)",
     "20301706"),
    ("chr7", None, None, "7p12 (GRB10), 7q32 (MEST)", "Silver-Russell syndrome", "-", "27585961"),
    ("chr11", 0, 22000000, "11p15 (H19/IGF2, KCNQ1OT1)", "Silver-Russell syndrome", "Beckwith-Wiedemann syndrome",
     "29377879;27585961"),
    ("chr14", None, None, "14q32.2 (DLK1-MEG3)", "Temple syndrome", "Kagami-Ogata syndrome", "24891339;26377239"),
    ("chr15", 20500000, 33400000, "15q11-q13 (SNRPN, UBE3A)", "Prader-Willi syndrome", "Angelman syndrome",
     "22237428;27615419"),
    ("chr20", None, None, "20q13.32 (GNAS)", "Mulchandani-Bhoj-Conlin syndrome (growth failure)",
     "Pseudohypoparathyroidism type 1B", "26248010;29959430"),
]


def autosome_unmasked_bp(bands=None):
    """AUTOSOME_BP without the autosomal MASKED_BANDS."""
    bands = masked_bands() if bands is None else bands
    return AUTOSOME_BP - sum(b - a for c, spans in bands.items() if c[3:] in AUTOSOMES for a, b in spans)


def roh_band(f):
    """The row of docs/11-roh-analysis.md's table an F_ROH is nearest to."""
    return next(label for top, label in ROH_BANDS if f < top)


def imprinted_rows(chrom, segs):
    """The IMPRINTED_REGIONS rows a flagged chromosome joins: its own
    chromosome's rows, and of a row with a region only when one of the runs
    (1-based, inclusive start-end) overlaps it."""
    out = []
    for c, a, b, region, mat, pat, pmid in IMPRINTED_REGIONS:
        if c != chrom:
            continue
        if a is not None and not any(masked_overlap({c: [(a, b)]}, c, s, e) for s, e in segs):
            continue
        out.append({"region": region, "maternal_upd": mat, "paternal_upd": pat, "pmid": pmid})
    return out


def sec_roh(d, s):
    """Step 11's bcftools roh output. Segments on chrX and chrY are left out
    (a man's single X is one long run), segments on any other contig that is
    not chr1 to chr22 (chrM, chrEBV, chrUn_*) are not counted, and so is the
    part of a segment that lies in a centromere or heterochromatin band
    (MASKED_BANDS), where few markers map. Total, Largest and the list of segments over 5 Mb are of the
    autosomal length that is left."""
    p = first_existing(d, [f"vcf/{s}_roh.txt", f"roh/{s}_roh.txt", "roh/*_roh.txt"])
    if not p:
        return None, {}
    bands = masked_bands()
    total = largest = masked = 0.0
    segs, sex_segs, big = 0, 0, []
    per_chrom = {}   # chrom -> [long bp, [(start, end), ...]] of the segments over ROH_MIN_BP
    with open(p, errors="replace") as f:
        for line in f:
            c = line.split()
            if len(c) < 6 or c[0] != "RG":
                continue
            try:
                start, end, length = int(c[3]), int(c[4]), float(c[5])
            except ValueError:
                continue
            if not math.isfinite(length):
                raise Unreadable(f"{os.path.basename(p)} has a run of length {c[5]}")
            if c[2] in ("chrX", "chrY", "X", "Y"):
                sex_segs += 1
                continue
            if not (c[2].startswith("chr") and c[2][3:] in AUTOSOMES):
                continue
            inside = min(length, masked_overlap(bands, c[2], start, end))
            masked += inside
            length -= inside
            if length <= 0:
                continue
            segs += 1
            total += length
            largest = max(largest, length)
            if length >= ROH_MIN_BP:
                pc = per_chrom.setdefault(c[2], [0.0, []])
                pc[0] += length
                pc[1].append((start, end))
                row = {"region": f"{c[2]}:{c[3]}-{c[4]}", "mb": round(length / 1e6, 1)}
                if inside:
                    row["masked_mb"] = round(inside / 1e6, 1)
                big.append(row)
    long_bp = sum(v[0] for v in per_chrom.values())
    f_roh = long_bp / autosome_unmasked_bp(bands)
    upd = None
    if per_chrom:
        chrom, (bp, spans) = max(per_chrom.items(), key=lambda kv: kv[1][0])
        if bp >= UPD_MIN_BP and bp > UPD_SHARE * long_bp:
            upd = {"chrom": chrom, "mb": round(bp / 1e6, 1), "share": round(bp / long_bp, 2),
                   "imprinted": imprinted_rows(chrom, spans)}
    return p, {"segments": segs, "total_mb": round(total / 1e6, 1), "largest_mb": round(largest / 1e6, 1),
               "autosomal_over_5mb": big, "masked_mb": round(masked / 1e6, 1), "sex_chromosome_segments": sex_segs,
               "long_mb": round(long_bp / 1e6, 1), "f_roh": round(f_roh, 4),
               # a chromosome that may come from one parent fills F_ROH by itself, so no band then
               "relationship": None if upd else roh_band(f_roh), "upd_flag": upd}


# Pages that show more of a called branch, free and with no upload. The URL
# patterns were checked live on 2026-10-25 (docs/37-y-haplogroup.md,
# "Explore your branch further"): FamilyTreeDNA Discover answers its own
# HaplogroupSearch for R-M269 and H5a7, and YFull's tree pages answer 200 for
# a branch it has and 404 for one it does not.
EXPLORE_DOCS = "https://github.com/GeiserX/OwnGenome/blob/main/docs/"


def haplogroup_explore(kind, hg):
    """[{name, url}] for a called branch: kind 'y' (Yleaf, e.g. R-M269 or
    R-M269*(xR-L23)) or 'mt' (haplogrep3, e.g. H5a7). Empty when there is no call."""
    from urllib.parse import quote
    hg = (hg or "").strip()
    if kind == "y":
        hg = hg.split("*", 1)[0]   # Yleaf's paragroup mark: the branch is the name before '*'
    if not hg or hg in (".", "NA", "insufficient markers"):
        return []
    q = quote(hg, safe="")
    if kind == "y":
        return [{"name": "FamilyTreeDNA Discover (free, no upload)", "url": f"https://discover.familytreedna.com/y-dna/{q}"},
                {"name": "YFull YTree", "url": f"https://www.yfull.com/tree/{q}/"},
                {"name": "More routes, and the files to share", "url": EXPLORE_DOCS + "37-y-haplogroup.md#explore-your-branch-further"}]
    return [{"name": "FamilyTreeDNA Discover (free, no upload)", "url": f"https://discover.familytreedna.com/mtdna/{q}"},
            {"name": "YFull MTree", "url": f"https://www.yfull.com/mtree/{q}/"},
            {"name": "More routes, and the files to share", "url": EXPLORE_DOCS + "12-mito-haplogroup.md#explore-your-branch-further"}]


def sec_haplogroup(d, s):
    p = first_existing(d, [f"mito/{s}_haplogroup.txt"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    r = rows[0] if rows else {}
    hg = (r.get("Haplogroup") or (list(r.values())[1] if len(r) > 1 else "") or "").strip('"')
    out = {"haplogroup": hg or ".", "explore": haplogroup_explore("mt", hg)}
    # haplocheck (step 12 with step 20's Mutect2 calls): a second haplogroup
    # in the allele fractions is another person's DNA in the sample.
    hc = first_existing(d, [f"mito/{s}_haplocheck.txt"])
    if hc:
        rows = [{k.strip('"'): (v or "").strip('"') for k, v in row.items() if k} for row in read_tsv(hc)]
        if not rows or "Contamination Status" not in rows[0]:
            raise Unreadable(f"{os.path.basename(hc)} has no Contamination Status column")
        out["contamination_status"] = rows[0]["Contamination Status"] or "."
        out["contamination_level"] = rows[0].get("Contamination Level") or "."
    return p, out


def sec_y_haplogroup(d, s):
    p = first_existing(d, [f"y_haplogroup/{s}_y_haplogroup.txt"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    if not rows or "Hg" not in rows[0]:
        raise Unreadable(f"{os.path.basename(p)} has no Hg column or no sample row")
    r = rows[0]
    hg = (r.get("Hg") or "").strip()
    return p, {"haplogroup": hg if hg not in ("", "NA") else "insufficient markers",
               "valid_markers": r.get("Valid_markers") or ".", "qc_score": r.get("QC-score") or ".",
               "explore": haplogroup_explore("y", hg)}


def sec_y_tree(d, s):
    """Step 37's placement on the current public tree (y_haplogroup/<id>_y_tree.txt,
    bin/y_tree_place.py): the terminal branch, how deep it sits, the SNP
    calls along the path, and the tree release the names come from. The
    per-level table beside it gives each branch of that path."""
    p = first_existing(d, [f"y_haplogroup/{s}_y_tree.txt"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    need = ("Hg", "Level", "Tree")
    if not rows or any(k not in rows[0] for k in need):
        raise Unreadable(f"{os.path.basename(p)} lacks a column of {', '.join(need)}")
    r = rows[0]
    hg = (r.get("Hg") or "").strip()
    out = {"haplogroup": hg if hg not in ("", "NA") else "no branch has a derived call",
           "level": r.get("Level") or ".", "tree": r.get("Tree") or ".",
           "path_derived": r.get("Path_derived") or ".", "path_ancestral": r.get("Path_ancestral") or ".",
           "path_no_call": r.get("Path_no_call") or ".",
           "ancestral_in_terminal": r.get("Ancestral_in_terminal") or ".",
           "outvoted_ancestral": r.get("Outvoted_ancestral") or ".",
           "snps_called": r.get("SNPs_called") or ".", "min_depth": r.get("Min_depth") or ".",
           "explore": haplogroup_explore("y", hg)}
    lvl = first_existing(d, [f"y_haplogroup/{s}_y_tree_path.tsv"])
    if lvl:
        # bin/y_tree_place.py writes the header alone when no branch is derived:
        # an empty path, not an unreadable file.
        with open_text(lvl) as f:
            rd = csv.DictReader(f, delimiter="\t")
            prows = list(rd)
            header = rd.fieldnames or []
        if "Branch" not in header:
            raise Unreadable(f"{os.path.basename(lvl)} has no Branch column")
        out["path"] = [{k: (row.get(k) or "") for k in
                        ("Level", "Branch", "SNPs", "On_tables", "Derived", "Ancestral", "No_call",
                         "Ancestral_SNPs")} for row in prows]
    return p, out


def sec_y_str(d, s):
    """Step 45's table (y_str/<id>_y_str.tsv, bin/ystr_report.py): one row per
    FamilyTreeDNA marker with its count and confidence, and the counts per
    confidence."""
    p = first_existing(d, [f"y_str/{s}_y_str.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    need = ("order", "panel", "marker", "count", "confidence")
    if not rows or any(k not in rows[0] for k in need):
        raise Unreadable(f"{os.path.basename(p)} has no rows or lacks a column of {', '.join(need)}")
    keep = [{k: r.get(k) or "" for k in need} for r in rows]
    counts = {c: sum(1 for r in keep if r["confidence"] == c) for c in ("high", "low", "no call", "not callable")}
    if sum(counts.values()) != len(keep):
        raise Unreadable(f"{os.path.basename(p)} has a confidence other than high, low, no call or not callable")
    return p, {"markers": keep, **{k.replace(" ", "_"): v for k, v in counts.items()}}


def sec_consumer_export(d, s):
    """Step 40's summary table (key, value): how many chip sites got a genotype."""
    p = first_existing(d, [f"consumer_export/{s}_consumer_export.tsv"])
    if not p:
        return None, {}
    data = {r["key"]: r.get("value", "") for r in read_tsv(p) if r.get("key")}
    for k in ("sites", "called", "no_call"):
        if not (data.get(k) or "").isdigit():
            raise Unreadable(f"{os.path.basename(p)} has no whole-number '{k}' row")
    reasons = {}
    for k, v in data.items():
        if k.startswith("no_call: "):
            if not (v or "").isdigit():
                raise Unreadable(f"{os.path.basename(p)} has a '{k}' row without a whole-number count")
            reasons[k[len("no_call: "):]] = int(v)
    return p, {"sites": int(data["sites"]), "called": int(data["called"]), "no_call": int(data["no_call"]),
               "from_reference_blocks": int(data.get("from_reference_blocks") or 0),
               "template_sites": data.get("template_sites") or ".",
               "file_23andme": data.get("file_23andme") or ".",
               "file_ancestrydna": data.get("file_ancestrydna") or ".",
               "no_call_reasons": reasons}


def sec_bloodtype(d, s):
    """Step 41's table: one row per blood group system, predicted from the genome."""
    p = first_existing(d, [f"bloodtype/{s}_bloodtype.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    if not rows or not {"system", "phenotype", "confidence"} <= set(rows[0]):
        raise Unreadable(f"{os.path.basename(p)} has no system, phenotype and confidence columns, or no row")
    systems = [{"system": r["system"], "phenotype": r["phenotype"] or ".", "genotype": r.get("genotype") or ".",
                "confidence": r["confidence"] or "."} for r in rows]
    by = {x["system"]: x["phenotype"] for x in systems}
    return p, {"systems": systems, "abo": by.get("ABO", "."), "rhd": by.get("RhD", ".")}


def sec_archaic(d, s):
    """Step 43's summary (archaic/<id>_archaic.tsv, bin/archaic_report.py): the
    Neanderthal and Denisovan shares and segment counts, the analysed span and
    the ten longest segments."""
    p = first_existing(d, [f"archaic/{s}_archaic.tsv"])
    if not p:
        return None, {}
    data = {r["key"]: r.get("value") or "" for r in read_tsv(p) if r.get("key")}  # a short row: None
    out = {"population": data.get("population", "."),
           "longest": [data[k] for k in sorted((k for k in data if re.fullmatch(r"longest_[0-9]+", k)),
                                               key=lambda k: int(k.split("_")[1]))]}
    for k in ("neanderthal_segments", "neanderthal_bp", "denisovan_segments", "denisovan_bp", "analysed_bp"):
        if not re.fullmatch(r"[0-9]+", data.get(k, "")):
            raise Unreadable(f"{os.path.basename(p)} has no whole number for {k}")
        out[k] = int(data[k])
    for k in ("neanderthal_percent", "denisovan_percent"):
        if not re.fullmatch(r"[0-9]+\.[0-9]+", data.get(k, "")):
            raise Unreadable(f"{os.path.basename(p)} has no percentage for {k}")
        out[k] = data[k]
    return p, out


# The two LPA risk variants of assets/traits.tsv (rsid -> its alt, the allele
# that goes with a high Lp(a)); tests/test_lpa_kiv2.py checks they agree.
LPA_RISK_ALLELES = {"rs10455872": "G", "rs3798220": "C"}


def sec_lpa(d, s):
    """Step 46's KIV-2 count (lpa/<id>_lpa_kiv2.tsv, bin/lpa_kiv2.py), with
    step 42's two LPA risk variants when the traits table is there, and
    whether either points to a genetically high Lp(a)."""
    p = first_existing(d, [f"lpa/{s}_lpa_kiv2.tsv"])
    if not p:
        return None, {}
    kv = {r.get("metric"): r.get("value") for r in read_tsv(p)}
    if kv.get("status") not in ("estimated", "too little depth", "not assessed"):
        raise Unreadable(f"{os.path.basename(p)} has no status row (estimated, too little depth or not assessed)")
    data = {k: kv.get(k) or "." for k in ("status", "kiv2_copies", "ci95_low", "ci95_high", "class",
                                          "kiv2_depth", "flank_depth", "mapq0_share", "message", "advice")}
    variants = []
    t = first_existing(d, [f"traits/{s}_traits.tsv"])
    for r in (read_tsv(t) if t else []):
        alt = LPA_RISK_ALLELES.get(r.get("rsid"))
        if not alt:
            continue
        called = r.get("status") == "called"
        n = re.split(r"[/|]", r.get("genotype") or "").count(alt) if called else None
        variants.append({"rsid": r["rsid"], "genotype": r.get("genotype") or ".", "status": r.get("status") or ".",
                         "risk_alleles": n})
    risk = sum(v["risk_alleles"] for v in variants if v["risk_alleles"] is not None)
    reasons = ([f"{risk} risk allele{'s' if risk != 1 else ''} at rs10455872 and rs3798220"] if risk else []) \
        + (["a low KIV-2 repeat count"] if data["class"] == "low" else [])
    read = data["status"] == "estimated" or any(v["risk_alleles"] is not None for v in variants)
    data.update(variants=variants, risk_alleles=risk, flag="high" if reasons else ".",
                flag_text=("Genetic signs of a high Lp(a): " + " and ".join(reasons) + ". Measure Lp(a) in blood."
                           if reasons else
                           "No genetic sign of a high Lp(a) from these markers. They explain only part of it: "
                           "the KIV-2 count explained 21% to 27% of the variation in a Danish population study."
                           if read else
                           "Not determined: neither the KIV-2 count nor the two risk variants could be read. "
                           "Measure Lp(a) in blood."))
    return p, data


VNTYPER_STATUSES = ("positive", "possible", "negative", "not assessed")


def sec_vntyper(d, s):
    """Step 52's table (vntyper/<id>_vntyper.tsv, bin/vntyper_report.py):
    VNtyper's MUC1 VNTR reading with its depth and sensitivity line."""
    p = first_existing(d, [f"vntyper/{s}_vntyper.tsv"])
    if not p:
        return None, {}
    kv = {r.get("metric"): r.get("value") for r in read_tsv(p)}
    if kv.get("status") not in VNTYPER_STATUSES:
        raise Unreadable(f"{os.path.basename(p)} has no status row ({', '.join(VNTYPER_STATUSES)})")
    return p, {k: kv.get(k) or "." for k in ("status", "confidence", "variant", "flag", "alt_depth",
                                             "vntr_mean_depth", "coverage_qc", "sensitivity", "message", "advice")}


LOCITYPER_COLUMNS = ("locus", "region", "status", "genotype", "haplotype_1", "haplotype_2", "quality",
                     "total_reads", "unexplained_reads", "weighted_distance", "warnings", "filter", "reason")


def sec_locityper(d, s):
    """Step 53's table (locityper/<id>_locityper.tsv, bin/locityper_report.py):
    one row per locus of assets/locityper_loci.bed, the two pangenome
    haplotypes Locityper picked with its quality and the authors' filter.
    collect() adds each locus's primary call (locityper_primary)."""
    p = first_existing(d, [f"locityper/{s}_locityper.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    if not rows or not set(LOCITYPER_COLUMNS) <= set(rows[0]):
        raise Unreadable(f"{os.path.basename(p)} lacks a row or one of the columns {', '.join(LOCITYPER_COLUMNS)}")
    loci = [{k: (r.get(k) or ".") for k in LOCITYPER_COLUMNS} for r in rows]
    return p, {"loci": loci, "called": sum(x["status"] == "called" for x in loci),
               "passed": sum(x["filter"] == "pass" for x in loci)}


def locityper_primary(sections):
    """The primary caller's call at each Locityper locus, as text, so the
    report shows the two side by side: T1K (step 08) for HLA-DRB1 and
    HLA-DQB1, T1K's KIR pass (KIR=true) for KIR, pypgx (step 32) and Cyrius
    (step 21) for CYP2D6. FCGR has no primary caller here. Locityper names
    pangenome haplotypes, not alleles, so the two are shown, not compared."""
    def old(key):  # a table older than the last run-all.sh run says so
        return " (stale: from an earlier run)" if sections[key]["state"] == "stale" else ""

    hla = sections["hla"]["data"] if sections["hla"]["state"] in ("ok", "stale") else {}
    by_gene = {x["gene"]: x for x in hla.get("loci", [])}
    out = {}
    for g in ("HLA-DRB1", "HLA-DQB1"):
        x = by_gene.get(g) or by_gene.get(g.replace("HLA-", ""))
        out[g] = ("T1K (step 08): " + ("/".join(x["alleles"]) or "no allele")
                  + (" (low confidence)" if x.get("low_confidence") else "") + old("hla")) if x else "T1K (step 08) did not run"
    kir = hla.get("kir")
    out["KIR"] = (f"T1K KIR pass (step 08): {len(kir['typed'])} genes typed{old('hla')}" if kir and kir.get("typed")
                  else "T1K's KIR pass (step 08, KIR=true) gave no call")
    cyp = []
    if sections["pypgx"]["state"] in ("ok", "stale"):
        cyp.append(f"pypgx (step 32): {sections['pypgx']['data'].get('cyp2d6', 'not called')}{old('pypgx')}")
    if sections["cyrius"]["state"] in ("ok", "stale"):
        cy = sections["cyrius"]["data"]
        usable = "" if cy.get("filter") == "PASS" else f" (not usable: Filter {cy.get('filter') or 'missing'})"
        cyp.append(f"Cyrius (step 21): {cy.get('genotype', 'none')}{usable}{old('cyrius')}")
    out["CYP2D6"] = "; ".join(cyp) or "pypgx (step 32) and Cyrius (step 21) did not run"
    out["FCGR"] = "no other step calls FCGR"
    return out


def sec_viral(d, s):
    """Step 47's table (viral/<id>_viral.tsv, bin/viral_report.py): one row per
    viral genome; the HHV-6 line and the genomes with reads (phiX included)."""
    p = first_existing(d, [f"viral/{s}_viral.tsv"])
    if not p:
        return None, {}
    rows = read_tsv(p)
    need = ("virus", "accession", "reads", "breadth", "relative_depth", "copies_per_cell", "call")
    if not rows or any(k not in rows[0] for k in need):
        raise Unreadable(f"{os.path.basename(p)} has no rows or lacks a column of {', '.join(need)}")
    keep = [{k: r[k] or "." for k in need + ("human_depth", "note")} for r in rows]
    ici = [r["virus"] for r in keep if r["call"] == "inherited integrated HHV-6 likely"]
    # No read on chr1 to chr22 leaves no human depth, so integration was never assessed.
    hhv6 = (f"inherited integrated {ici[0]} likely" if ici
            else "not assessed (no human depth)" if keep[0]["human_depth"] == "."
            else "no integrated HHV-6")
    return p, {"hhv6": hhv6,
               "human_depth": keep[0]["human_depth"],
               "found": [r for r in keep if r["reads"].isdigit() and int(r["reads"]) > 0], "genomes": len(keep)}


def sec_chip(d, s):
    """Step 49's tables (chip/<id>_chip.tsv and chip/<id>_chip_depth.tsv, bin/chip_report.py):
    the CHIP calls, how many other protein changes were listed and not counted,
    and the depth and detection limit per gene with their median over the
    covered genes. The limit is the one at the median depth, as step 49 prints it."""
    p = first_existing(d, [f"chip/{s}_chip.tsv"])
    dp = first_existing(d, [f"chip/{s}_chip_depth.tsv"])
    if not p or not dp:
        return None, {}
    calls, depth = read_tsv(p), read_tsv(dp)
    need = ("gene", "protein_change", "vaf", "alt_reads", "depth", "status", "note")
    with open(p, encoding="utf-8", errors="replace") as f:
        head = f.readline().rstrip("\n").split("\t")
    if any(k not in head for k in need):
        raise Unreadable(f"{os.path.basename(p)} lacks a column of {', '.join(need)}")
    if not depth or not {"gene", "mean_depth", "detection_limit_vaf", "status"} <= set(depth[0]):
        raise Unreadable(f"{os.path.basename(dp)} has no gene, mean_depth, detection_limit_vaf and status columns, or no row")
    chip = [{k: r.get(k) or "." for k in need} for r in calls if r.get("status") == "chip"]
    covered = sorted(float(r["mean_depth"]) for r in depth if r["status"] == "covered")
    med = None
    if covered:
        n = len(covered)
        med = covered[n // 2] if n % 2 else (covered[n // 2 - 1] + covered[n // 2]) / 2
    return p, {"calls": chip, "other": len(calls) - len(chip),
               "genes": [{k: r.get(k) or "." for k in ("gene", "mean_depth", "detection_limit_vaf", "status")}
                         for r in depth],
               "median_depth": med, "median_limit": chip_report.detection_limit(med) if med is not None else None,
               "not_covered": [r["gene"] for r in depth if r["status"] != "covered"]}


def sec_mei(d, s):
    """Step 50's tables (mei/<id>_mei.tsv and _mei_sites.tsv, bin/mei_report.py):
    the insertions MEGAnE calls, counted per element, and every row of the
    known-site catalog with seen, not seen or not callable."""
    p = first_existing(d, [f"mei/{s}_mei_sites.tsv"])
    if not p:
        return None, {}
    sites = read_tsv(p)
    need = ("gene", "element", "condition", "inheritance", "status", "mean_depth", "call")
    if not sites or any(k not in sites[0] for k in need):
        raise Unreadable(f"{os.path.basename(p)} has no rows or lacks a column of {', '.join(need)}")
    cp = os.path.join(os.path.dirname(p), f"{s}_mei.tsv")
    if not os.path.isfile(cp):
        raise Unreadable(f"{os.path.basename(cp)} is missing beside {os.path.basename(p)}")
    with open_text(cp) as f:   # the header alone: a run with no insertion writes no row
        r = csv.DictReader(f, delimiter="\t")
        head, calls = r.fieldnames or [], list(r)
    if any(k not in head for k in ("chrom", "pos", "element", "copies")):
        raise Unreadable(f"{os.path.basename(cp)} lacks a column of chrom, pos, element, copies")
    counts = {}
    for c in calls:
        counts[c["element"]] = counts.get(c["element"], 0) + 1
    return p, {"calls": len(calls), "by_element": dict(sorted(counts.items())),
               "sites": [{k: r[k] or "." for k in need} for r in sites]}


MOCHA_SUMMARY_KEYS = ("status", "message", "depth", "het_sites", "het_per_10mb", "chromosomes", "floor_read",
                      "floor_min", "floor", "calls_mocha", "calls_reported", "phasing")
MOCHA_CALL_KEYS = ("chrom", "start", "end", "length_mb", "arms", "type", "cell_fraction", "n_hets", "lod_baf_phase",
                   "reported", "reason")


def sec_mocha(d, s):
    """Step 51's tables (mosaic/<id>_mocha_summary.tsv and _mocha.tsv,
    bin/mocha_report.py): the status, the floor, and the calls reported."""
    p = first_existing(d, [f"mosaic/{s}_mocha_summary.tsv"])
    if not p:
        return None, {}
    kv = {r.get("metric"): r.get("value") for r in read_tsv(p)}
    if kv.get("status") not in ("assessed", "not assessed"):
        raise Unreadable(f"{os.path.basename(p)} has no status row (assessed or not assessed)")
    data = {k: kv.get(k) or "." for k in MOCHA_SUMMARY_KEYS}
    c = first_existing(d, [f"mosaic/{s}_mocha.tsv"])
    if not c:
        raise Unreadable(f"{os.path.basename(p)} has no {s}_mocha.tsv beside it")
    with open_text(c) as f:
        head = f.readline().rstrip("\r\n").split("\t")
    if any(k not in head for k in MOCHA_CALL_KEYS):
        raise Unreadable(f"{os.path.basename(c)} lacks a column of {', '.join(MOCHA_CALL_KEYS)}")
    rows = read_tsv(c)
    data["calls"] = [{k: r[k] or "." for k in MOCHA_CALL_KEYS} for r in rows if r.get("reported") == "yes"]
    data["not_reported"] = sum(1 for r in rows if r.get("reported") != "yes")
    if data["calls_reported"] != "." and data["calls_reported"] != str(len(data["calls"])):
        raise Unreadable(f"{os.path.basename(c)} has {len(data['calls'])} reported calls, "
                         f"{os.path.basename(p)} says {data['calls_reported']}")
    return p, data


def sec_traits(d, s):
    """Step 42's traits table (traits/<id>_traits.tsv, bin/traits_report.py):
    every row with its genotype, status and reading, the counts per status,
    and the APOE line on its own (the reports fold it away)."""
    p = first_existing(d, [f"traits/{s}_traits.tsv"])
    if not p:
        return None, {}
    with open(p, encoding="utf-8", errors="replace") as f:
        rows = list(csv.DictReader(f, delimiter="\t", quoting=csv.QUOTE_NONE))
    need = ("rsid", "gene", "trait", "genotype", "status", "reading", "ancestry", "note", "pmids")
    if not rows or any(k not in rows[0] for k in need):
        raise Unreadable(f"{os.path.basename(p)} has no rows or lacks a column of {', '.join(need)}")
    # A row cut short (DictReader fills the missing fields with None) or with
    # extra fields (kept under the key None) is a damaged table, never a
    # called trait without its reading.
    bad = next((i for i, r in enumerate(rows, 2) if None in r or any(r.get(k) is None for k in need)), None)
    if bad:
        raise Unreadable(f"{os.path.basename(p)} line {bad} does not have one field per column")
    keep = [{k: r[k] for k in need} for r in rows]
    apoe = [r for r in keep if r["gene"] == "APOE"]
    counts = {st: sum(1 for r in keep if r["status"] == st and r["gene"] != "APOE")
              for st in ("called", "not_covered", "other_allele")}
    return p, {"rows": [r for r in keep if r["gene"] != "APOE"], "apoe": apoe, **counts,
               "rules": genotype_rules(d, s)}


def genotype_rules(d, s):
    """Step 42's genotype rules (traits/<id>_genotype_rules.tsv): one line per
    rule that reads several sites together (TAS2R38's haplotypes), or None
    when the table is not there (a run from before the rules). The rules
    whose line the traits table already holds (APOE) are left out, so the
    reports keep folding APOE away."""
    p = first_existing(d, [f"traits/{s}_genotype_rules.tsv"])
    if not p:
        return None
    with open(p, encoding="utf-8", errors="replace") as f:
        rows = list(csv.DictReader(f, delimiter="\t", quoting=csv.QUOTE_NONE))
    need = traits_report.RULE_OUT_COLUMNS
    if rows and any(k not in rows[0] for k in need):
        raise Unreadable(f"{os.path.basename(p)} lacks a column of {', '.join(need)}")
    bad = next((i for i, r in enumerate(rows, 2) if None in r or any(r.get(k) is None for k in need)), None)
    if bad:
        raise Unreadable(f"{os.path.basename(p)} line {bad} does not have one field per column")
    out = [{k: r[k] for k in need} for r in rows if r["rule"] not in traits_report.TRAITS_TABLE_RULES]
    for r in out:
        if r["rule"] == genotype_rules_indels.G2["rule"]:
            g2 = apol1_g2(d, s)
            r["g2"] = dict(g2, **genotype_rules_indels.combine(r, g2))
    return out


def vcf_window(path, chrom, start, end):
    """The split data lines of a VCF on chrom that can cover start..end: read
    from the tabix linear index's offset for start when the .tbi is there,
    else the whole file (stopping past the window)."""
    out = []
    if os.path.exists(path + ".tbi"):
        ioff = tabix_linear_index(path).get(chrom)
        if not ioff:
            return out
        voff = ioff[min((start - 1) >> 14, len(ioff) - 1)]
        voff = 0 if voff == 0xFFFFFFFFFFFFFFFF else voff
        with open(path, "rb") as raw:
            raw.seek(voff >> 16)
            g = gzip.GzipFile(fileobj=raw)
            g.read(voff & 0xFFFF)
            lines = (x.decode(errors="replace").rstrip("\n").split("\t") for x in g if not x.startswith(b"#"))
            seen = False
            for r in lines:
                if len(r) < 2 or not r[1].isdigit():
                    continue
                if r[0] != chrom:
                    if seen:
                        break
                    continue
                seen = True
                if int(r[1]) > end:
                    break
                out.append(r)
        return out
    seen = False
    for r in vcf_records(path):
        if len(r) < 2 or not r[1].isdigit():
            continue
        if r[0] != chrom:
            if seen:
                break
            continue
        seen = True
        if int(r[1]) > end:
            break
        out.append(r)
    return out


def apol1_g2(d, s):
    """APOL1 G2 (rs71785313, a 6-base deletion step 42's rule table cannot
    read) from the sample's gVCF when it is there (vcf/<id>.g.vcf.gz, whose
    reference blocks tell a reference base from no reads), else its VCF
    (vcf/<id>.vcf.gz): bin/genotype_rules_indels.py's read_g2, plus source."""
    m = genotype_rules_indels
    for name, kind in ((f"vcf/{s}.g.vcf.gz", "gvcf"), (f"vcf/{s}.vcf.gz", "vcf")):
        p = os.path.join(d, name)
        if not os.path.isfile(p):
            continue
        try:
            recs = vcf_window(p, m.G2["chrom"], m.CONTEXT_START, m.SPAN[1])
        except (OSError, EOFError, ValueError, struct.error, gzip.BadGzipFile, zlib.error, Unreadable) as e:
            return {"status": "not_read", "copies": None, "gq": ".", "source": name,
                    "text": f"G2 (rs71785313) not read: {name} could not be read ({e})."}
        return dict(m.read_g2(recs, kind), source=name)
    return {"status": "not_read", "copies": None, "gq": ".", "source": None,
            "text": "G2 (rs71785313) not read: no VCF of the sample in this folder."}


def sec_carrier(d, s):
    """The sample's own carrier list (carrier_list) from step 06's ClinVar
    hits; its own section, so the ClinVar section reads the same whether or
    not the assets folder is there."""
    p = first_existing(d, [f"clinvar/{s}_clinvar_hits.vcf"])
    if not p:
        return None, {}
    return p, carrier_list_or_note(p, d, s)


def sec_exomiser(d, s):
    """Step 48's ranking (exomiser/<id>_exomiser_top.tsv, bin/exomiser_report.py):
    the HPO terms searched, and each rank with its gene, mode of inheritance,
    scores and contributing variants. A table without rows means no gene
    passed Exomiser's filters."""
    p = first_existing(d, [f"exomiser/{s}_exomiser_top.tsv"])
    if not p:
        return None, {}
    with open(p, encoding="utf-8", errors="replace") as f:
        rows = list(csv.DictReader(f, delimiter="\t", quoting=csv.QUOTE_NONE))
        f.seek(0)
        header = f.readline().rstrip("\n").split("\t")
    need = ("rank", "gene", "moi", "exomiser_score", "phenotype_score", "variant", "genotype", "consequence",
            "hgvs", "acmg", "disease")
    if any(k not in header for k in need):
        raise Unreadable(f"{os.path.basename(p)} lacks a column of {', '.join(need)}")
    if any(None in r or any(r.get(k) is None for k in need) for r in rows):
        raise Unreadable(f"{os.path.basename(p)} has a line that does not have one field per column")
    terms = []
    tp = os.path.join(os.path.dirname(p), f"{s}_exomiser_terms.tsv")
    if os.path.isfile(tp):
        with open(tp, encoding="utf-8", errors="replace") as f:
            terms = [{"id": r.get("hpo_id", ""), "label": r.get("label", "")}
                     for r in csv.DictReader(f, delimiter="\t", quoting=csv.QUOTE_NONE)]
    genes = []
    for r in rows:
        if not genes or genes[-1]["rank"] != r["rank"]:
            genes.append({k: r[k] for k in ("rank", "gene", "moi", "exomiser_score", "phenotype_score")})
            genes[-1]["variants"] = []
        genes[-1]["variants"].append({k: r[k] for k in ("variant", "genotype", "consequence", "hgvs", "acmg", "disease")})
    return p, {"terms": terms, "genes": genes}


def sec_couple(d, s):
    """Step 39's couple table (couple/<id>_couple_carrier.tsv, bin/couple_report.py):
    the partner, the shared and X-linked genes with the chance for a child, and
    how many genes only one partner carries a change in. A table without rows
    (no finding in either partner) names no partner."""
    p = first_existing(d, [f"couple/{s}_couple_carrier.tsv"])
    if not p:
        return None, {}
    with open(p, encoding="utf-8", errors="replace") as f:
        reader = csv.DictReader(f, delimiter="\t", quoting=csv.QUOTE_NONE)
        rows = list(reader)
        cols = reader.fieldnames or []
    need = ("gene", "condition", "result", "offspring_risk", "partner_a", "findings_a", "partner_b",
            "findings_b", "explanation")
    if any(k not in cols for k in need):
        raise Unreadable(f"{os.path.basename(p)} lacks a column of {', '.join(need)}")
    bad = next((i for i, r in enumerate(rows, 2) if None in r or any(r.get(k) is None for k in need)
                or s not in (r["partner_a"], r["partner_b"]) or r["result"] not in ("shared", "x_linked", "one_partner")),
               None)
    if bad:
        raise Unreadable(f"{os.path.basename(p)} line {bad} is not a row of {s}'s couple table")

    def side(r):
        own = "a" if r["partner_a"] == s else "b"
        other = "b" if own == "a" else "a"
        return {"gene": r["gene"], "condition": r["condition"], "risk": r["offspring_risk"],
                "own": r["findings_" + own], "partner_findings": r["findings_" + other],
                "explanation": r["explanation"]}
    partner = next((r["partner_b"] if r["partner_a"] == s else r["partner_a"] for r in rows), None)
    return p, {"partner": partner,
               "shared": [side(r) for r in rows if r["result"] == "shared"],
               "x_linked": [side(r) for r in rows if r["result"] == "x_linked"],
               "one_partner": sum(1 for r in rows if r["result"] == "one_partner")}


# The D-loop poly-C tracts (rCRS m.303-315 and m.16184-16193). Length
# variants there, an extra or missing C, are found in most people, often as a
# mix of lengths that reads as heteroplasmy; they are counted apart and left
# out of the heteroplasmic count.
POLYC_TRACTS = ((303, 315), (16184, 16193))


def mito_span(pos, ref, alt):
    """The rCRS positions a chrM allele changes: an SNV or MNV its own bases,
    an indel the bases after its anchor base (Mutect2 writes the anchor)."""
    if len(ref) == len(alt):
        return pos, pos + len(ref) - 1
    return pos + 1, pos + max(len(ref) - 1, 1)


def in_polyc(pos, ref, alt):
    lo, hi = mito_span(pos, ref, alt)
    return any(lo <= b and hi >= a for a, b in POLYC_TRACTS)


def min_allele(pos, ref, alt):
    """(pos, ref, alt) without the bases ref and alt share at the end and then
    at the start, one base always kept: Mutect2's A>AC of a two-ALT record
    AC>A,ACC becomes the same key as ClinVar's and gnomAD's."""
    while len(ref) > 1 and len(alt) > 1 and ref[-1] == alt[-1]:
        ref, alt = ref[:-1], alt[:-1]
    while len(ref) > 1 and len(alt) > 1 and ref[0] == alt[0]:
        ref, alt, pos = ref[1:], alt[1:], pos + 1
    return pos, ref, alt


MITO_ANNOT_COLS = ["chrom", "pos", "ref", "alt", "filter", "af", "clinvar_id", "clinvar_significance",
                   "clinvar_review", "clinvar_condition", "gnomad_af_hom", "gnomad_af_het", "gnomad_max_hl"]
# ClinVar classes: P/LP are counted; a record whose every class is benign or
# carries no classification is not listed unless the call is heteroplasmic.
CLINVAR_PLP = ("Pathogenic", "Likely_pathogenic", "Pathogenic/Likely_pathogenic")
CLINVAR_BENIGN = ("Benign", "Likely_benign", "Benign/Likely_benign", "not_provided")


def clinvar_terms(sig):
    return [t for t in re.split(r"[|,]", sig or "") if t and t != "."]


def vcf_file_date(path):
    with open_text(path) as f:
        for line in f:
            if not line.startswith("#"):
                break
            if line.startswith("##fileDate="):
                return line.strip().split("=", 1)[1]
    return None


def mito_annotate(vcf, clinvar=None, gnomad=None):
    """Step 20's annotation table: each ALT of each record of the chrM VCF,
    with its ClinVar record (CLNSIG, CLNREVSTAT, CLNDN of the same allele) and
    gnomAD v3.1's chrM sites (AF_hom, AF_het, max_hl of the same allele).
    Returns (header lines, rows); a file not given or not there is named in
    the header and its columns are '.'."""
    head, cv, gn = [], {}, {}
    if clinvar and os.path.isfile(clinvar) and os.path.isfile(clinvar + ".tbi"):
        head.append(f"# clinvar\t{os.path.basename(clinvar)}\t{vcf_file_date(clinvar) or 'fileDate unknown'}")
    else:
        head.append(f"# clinvar\tnot annotated\t{'no ' + clinvar + ' with its .tbi' if clinvar else 'no ClinVar file given'}")
        clinvar = None
    if gnomad and os.path.isfile(gnomad):
        head.append(f"# gnomad\t{os.path.basename(gnomad)}\t-")
        for r in vcf_records(gnomad):
            if len(r) < 8:
                continue
            i = info_map(r[7])
            for alt in r[4].split(","):
                gn[min_allele(int(r[1]), r[3], alt)] = (i.get("AF_hom", "."), i.get("AF_het", "."), i.get("max_hl", "."))
    else:
        head.append(f"# gnomad\tnot annotated\t{'no ' + gnomad if gnomad else 'no gnomAD chrM file given'}")
    recs = list(vcf_records(vcf))
    if clinvar:
        index = tabix_linear_index(clinvar)
        chrom = next((c for c in ("chrM", "MT") if c in index), None)
        wanted = {int(r[1]) for r in recs if len(r) >= 5} | {min_allele(int(r[1]), r[3], a)[0] for r in recs if len(r) >= 5 for a in r[4].split(",")}
        for pos in sorted(wanted) if chrom else ():
            # ClinVar's own record can start a base earlier or later than
            # Mutect2's for the same indel; the alleles decide.
            for q in (x for x in (pos - 1, pos, pos + 1) if x >= 1):
                for c in tabix_lines(clinvar, index, chrom, q):
                    i = info_map(c[7]) if len(c) > 7 else {}
                    for alt in c[4].split(","):
                        cv[min_allele(int(c[1]), c[3], alt)] = (c[2], i.get("CLNSIG", "."), i.get("CLNREVSTAT", "."),
                                                                i.get("CLNDN", "."))
    rows = []
    for r in recs:
        if len(r) < 8:
            continue
        fmt = r[8].split(":") if len(r) > 9 else []
        val = r[9].split(":") if len(r) > 9 else []
        afs = val[fmt.index("AF")].split(",") if "AF" in fmt and fmt.index("AF") < len(val) else []
        for k, alt in enumerate(r[4].split(",")):
            key = min_allele(int(r[1]), r[3], alt)
            rows.append([r[0], str(key[0]), key[1], key[2], r[6], afs[k] if k < len(afs) else ".",
                         *cv.get(key, (".", ".", ".", ".")), *gn.get(key, (".", ".", "."))])
    return head, rows


def mito_annotate_main(argv):
    ap = argparse.ArgumentParser(prog="collect_summary.py mito-annotate",
                                 description="Write step 20's chrM annotation table (ClinVar and gnomAD).")
    ap.add_argument("--vcf", required=True, help="the filtered chrM VCF (<sample>_chrM_filtered.vcf.gz)")
    ap.add_argument("--clinvar", default="", help="ClinVar's VCF with chr names, and its .tbi")
    ap.add_argument("--gnomad", default="", help="gnomAD v3.1's chrM sites VCF")
    ap.add_argument("--out", required=True)
    a = ap.parse_args(argv)
    try:
        head, rows = mito_annotate(a.vcf, a.clinvar or None, a.gnomad or None)
    except (Unreadable, OSError, ValueError, struct.error, zlib.error) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        return 1
    tmp = a.out + ".tmp"
    with open(tmp, "w") as f:
        f.write("".join(h + "\n" for h in head))
        f.write("\t".join(MITO_ANNOT_COLS) + "\n")
        f.write("".join("\t".join(r) + "\n" for r in rows))
    os.replace(tmp, a.out)
    for h in head:
        c = h[2:].split("\t")
        print(f"  {c[0]}: {c[1]}" + (f" ({c[2]})" if c[2] != "-" else ""))
    n_cv = sum(r[7] != "." for r in rows)
    n_gn = sum(r[10] != "." for r in rows)
    print(f"  {len(rows)} chrM alleles: {n_cv} with a ClinVar record, {n_gn} in gnomAD")
    return 0


# --- step 20: large mtDNA deletions (eKLIPse) and the mtDNA copy number ----------

# <id>_chrM_deletions.tsv: eKLIPse's eKLIPse_deletions.csv, one row per
# deletion. bkp5 and bkp3 are its 5' and 3' breakpoints as it writes them: the
# last rCRS base kept before the deletion and the first one kept after it, so
# the deleted bases are bkp5+1 to bkp3-1 (across the origin, m.16569 to m.1,
# when bkp3 is the smaller). level is the share of the mtDNA
# molecules with the deletion (in %), the mean of the forward and reverse
# strand estimates.
MITO_DEL_COLS = ["bkp5", "bkp3", "length", "level_pct", "level_fwd_pct", "level_rev_pct", "blast_5", "blast_3",
                 "depth_5", "depth_3", "repeat"]
# A deletion of assets/mito_confirmed.tsv is named when both breakpoints lie
# within this many bases of its own: a deletion between two direct repeats
# can be written anywhere along them (the common one's repeats are 13 bp).
MITO_DEL_TOLERANCE = 50
# The report lists every named deletion and any other at this level or more;
# the rest are counted.
MITO_DEL_LIST_PCT = 1.0
MITO_LENGTH = 16569  # the rCRS, GRCh38's chrM
# The mtDNA copy number: chrM depth over autosomal depth, from the read counts
# of the BAM index (samtools idxstats), times 2 for the two nuclear copies.
# Under this many reads per kb of autosome (about 3x at 150 bp) the BAM is a
# slice, a panel or a shallow genome, and the ratio would mean nothing. An
# exome can pass the floor and is not told apart from a genome here: the
# number is meant for whole-genome BAMs (docs/20-mtoolbox.md).
MITO_CN_MIN_READS_PER_KB = 20
MITO_CN_COLS = ["status", "reason", "chrm_reads", "chrm_length", "autosomal_reads", "autosomal_length",
                "depth_ratio", "mtdna_copies"]


def _eklipse_num(x):
    v = float(x.replace(",", "."))
    if not math.isfinite(v):
        raise ValueError(f"{x!r} is not a number")
    return v


def eklipse_deletions(path):
    """The rows of eKLIPse's eKLIPse_deletions.csv (';' between fields, each
    quoted, a comma for the decimal point), as MITO_DEL_COLS dicts."""
    with open(path, errors="replace", newline="") as f:
        rows = [r for r in csv.reader(f, delimiter=";", quotechar='"') if r]
    if not rows or rows[0][:3] != ["Title", "5' breakpoint", "3' breakpoint"]:
        raise ValueError(f"{path}: not eKLIPse's deletions table (header {rows[0][:3] if rows else 'missing'})")
    out = []
    for r in rows[1:]:
        if len(r) < 11:
            raise ValueError(f"{path}: a row with {len(r)} fields, want 11: {';'.join(r)}")
        b5, b3 = int(r[1]), int(r[2])
        out.append({"bkp5": b5, "bkp3": b3, "length": b3 - b5 - 1 if b3 > b5 else MITO_LENGTH - b5 + b3 - 1,
                    "level_pct": round(_eklipse_num(r[3]), 3), "level_fwd_pct": round(_eklipse_num(r[4]), 3),
                    "level_rev_pct": round(_eklipse_num(r[5]), 3), "blast_5": int(r[6]), "blast_3": int(r[7]),
                    "depth_5": int(_eklipse_num(r[8])), "depth_3": int(_eklipse_num(r[9])), "repeat": r[10] or "none"})
    return out


def mito_deletions_main(argv):
    ap = argparse.ArgumentParser(prog="collect_summary.py mito-deletions",
                                 description="Write step 20's large-deletion table from eKLIPse's deletions CSV.")
    ap.add_argument("--eklipse", required=True, help="eKLIPse's eKLIPse_deletions.csv")
    ap.add_argument("--out", required=True, help="<sample>_chrM_deletions.tsv")
    a = ap.parse_args(argv)
    try:
        rows = eklipse_deletions(a.eklipse)
    except (OSError, ValueError) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        return 1
    tmp = a.out + ".tmp"
    with open(tmp, "w") as f:
        f.write("\t".join(MITO_DEL_COLS) + "\n")
        for r in rows:
            f.write("\t".join("." if r[c] is None else str(r[c]) for c in MITO_DEL_COLS) + "\n")
    os.replace(tmp, a.out)
    print(f"  {len(rows)} large deletions (eKLIPse)")
    return 0


def mito_copies(idxstats):
    """The mtDNA copy number block from `samtools idxstats` text: chrM reads
    per base over autosomal reads per base, times 2."""
    chrm = None
    auto_reads = auto_len = 0
    for line in idxstats.splitlines():
        c = line.split("\t")
        if len(c) < 3 or c[0] == "*":
            continue
        name, length, mapped = c[0], int(c[1]), int(c[2])
        if name in ("chrM", "MT"):
            chrm = (mapped, length)
        elif re.fullmatch(r"(chr)?([1-9]|1[0-9]|2[0-2])", name):
            auto_reads += mapped
            auto_len += length
    b = dict.fromkeys(MITO_CN_COLS)
    b.update(status="not assessed", autosomal_reads=auto_reads, autosomal_length=auto_len)
    if chrm:
        b.update(chrm_reads=chrm[0], chrm_length=chrm[1])
    if not chrm or not chrm[0] or not chrm[1]:
        b["reason"] = "no_chrM_reads"
    elif not auto_len or auto_reads / auto_len * 1000 < MITO_CN_MIN_READS_PER_KB:
        b["reason"] = "low_autosomal_depth"
    else:
        ratio = (chrm[0] / chrm[1]) / (auto_reads / auto_len)
        b.update(status="called", depth_ratio=round(ratio, 3), mtdna_copies=round(2 * ratio, 1))
    return b


def mito_copies_main(argv):
    ap = argparse.ArgumentParser(prog="collect_summary.py mito-copies",
                                 description="Write step 20's mtDNA copy number from samtools idxstats.")
    ap.add_argument("--idxstats", required=True, help="samtools idxstats of the whole BAM")
    ap.add_argument("--out", required=True, help="<sample>_mtdna_copies.tsv")
    a = ap.parse_args(argv)
    try:
        with open(a.idxstats, errors="replace") as f:
            b = mito_copies(f.read())
    except (OSError, ValueError) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        return 1
    tmp = a.out + ".tmp"
    with open(tmp, "w") as f:
        f.write("metric\tvalue\n")
        f.write("".join(f"{k}\t{'.' if b[k] is None else b[k]}\n" for k in MITO_CN_COLS))
    os.replace(tmp, a.out)
    print(f"  mtDNA copy number: {b['mtdna_copies'] if b['status'] == 'called' else 'not assessed (' + b['reason'] + ')'}")
    return 0


def mito_confirmed():
    """(SNVs {(pos, ref, alt): row}, deletions [row]) of assets/mito_confirmed.tsv."""
    snv, dels = {}, []
    for r in asset_rows("mito_confirmed.tsv"):
        if r["kind"] == "snv":
            snv[(int(r["pos"]), r["ref"], r["alt"])] = r
        elif r["kind"] == "deletion":
            dels.append(r)
    return snv, dels


def known_deletion(row, dels):
    """The deletion of assets/mito_confirmed.tsv a called deletion is, or None."""
    for k in dels:
        if abs(row["bkp5"] + 1 - int(k["pos"])) <= MITO_DEL_TOLERANCE \
                and abs(row["bkp3"] - 1 - int(k["end"])) <= MITO_DEL_TOLERANCE:
            return k
    return None


def _confirmed_fields(k):
    return {"hgvs": k["hgvs"], "gene": k["gene"], "disease": k["disease"], "clinvar_id": k["clinvar_id"],
            "pmids": k["pmids"], "note": k["heteroplasmy_note"]}


def mito_extended(d, s, vcf):
    """Step 20's confirmed variants (assets/mito_confirmed.tsv) among the chrM
    calls, its large deletions and the mtDNA copy number."""
    out = {}
    try:
        snv, dels = mito_confirmed()
    except Unreadable as e:
        snv, dels = None, []
        out["confirmed_note"] = str(e)
    confirmed = []
    for r in vcf_records(vcf) if snv is not None else ():
        if len(r) < 8:
            continue
        fmt = r[8].split(":") if len(r) > 9 else []
        val = r[9].split(":") if len(r) > 9 else []
        afs = val[fmt.index("AF")].split(",") if "AF" in fmt and fmt.index("AF") < len(val) else []
        for i, alt in enumerate(r[4].split(",")):
            k = snv.get(min_allele(int(r[1]), r[3], alt))
            if k:
                confirmed.append(_confirmed_fields(k) | {"pos": int(k["pos"]), "filter": r[6],
                                                         "af": afs[i] if i < len(afs) else "."})
    if snv is not None:
        out["confirmed"] = confirmed
    p = first_existing(d, [f"mito/{s}_chrM_deletions.tsv"])
    if not p:
        out["deletions_note"] = ("no deletion table (mito/<id>_chrM_deletions.tsv): step 20 ran before it looked "
                                 "for large deletions, eKLIPse wrote no table (see the step's log), or the Nextflow "
                                 "MITO_DELETIONS task failed (see the Nextflow log)")
    else:
        with open(p, errors="replace") as f:
            rows = list(csv.DictReader(f, delimiter="\t"))
        listed = []
        for r in rows:
            row = {"bkp5": int(r["bkp5"]), "bkp3": int(r["bkp3"]),
                   "length": None if r["length"] == "." else int(r["length"]), "level_pct": float(r["level_pct"])}
            k = known_deletion(row, dels)
            if k or row["level_pct"] >= MITO_DEL_LIST_PCT:
                listed.append(row | ({"known": _confirmed_fields(k)} if k else {}))
        out["deletions"] = {"source": os.path.relpath(p, d), "total": len(rows), "listed": listed,
                            "list_floor_pct": MITO_DEL_LIST_PCT}
    p = first_existing(d, [f"mito/{s}_mtdna_copies.tsv"])
    if p:
        m = dict(read_tsv_pairs(p))
        out["copies"] = {k: (None if m.get(k) in (None, ".") else m[k]) for k in MITO_CN_COLS}
    else:
        out["copies_note"] = ("no copy-number table (mito/<id>_mtdna_copies.tsv): step 20 ran before it counted "
                              "mtDNA copies, or the Nextflow MITO_DELETIONS task failed (see the Nextflow log)")
    return out


def read_tsv_pairs(path):
    """[(metric, value)] of a two-column table with a header line."""
    with open(path, errors="replace") as f:
        lines = [ln.rstrip("\n").split("\t") for ln in f if ln.strip()]
    return [(c[0], c[1]) for c in lines[1:] if len(c) >= 2]


def mito_annotation(d, s, vcf):
    """Step 20's annotation table, as the report lists it, or a note on why
    there is none."""
    p = first_existing(d, [f"mito/{s}_chrM_annotated.tsv"])
    if not p:
        return {"note": "no annotation table (mito/<id>_chrM_annotated.tsv): step 20 ran before it annotated, "
                        "or the Nextflow MITO_VARIANTS_ANNOTATE task failed (see the Nextflow log)"}
    if os.path.getmtime(p) < os.path.getmtime(vcf):
        return {"note": f"{os.path.relpath(p, d)} is older than the chrM calls: rerun step 20"}
    sources, lines = {}, []
    with open(p, errors="replace") as f:
        for line in f:
            if line.startswith("# "):
                c = line[2:].rstrip("\n").split("\t")
                if len(c) >= 3:
                    sources[c[0]] = (f"not annotated: {c[2]}" if c[1] == "not annotated"
                                     else c[1] + (f" ({c[2]})" if c[2] != "-" else ""))
            elif line.strip():
                lines.append(line)
    rows = [r for r in csv.DictReader(lines, delimiter="\t") if r.get("filter") == "PASS"]
    return {"source": os.path.relpath(p, d), "clinvar": sources.get("clinvar"), "gnomad": sources.get("gnomad"),
            "rows": rows}


def sec_mito(d, s):
    p = first_existing(d, [f"mito/{s}_chrM_filtered.vcf.gz"])
    if not p:
        return None, {}
    passed = het = 0
    polyc, sites = [], []
    for r in vcf_records(p):
        if len(r) < 10 or r[6] != "PASS":
            continue
        passed += 1
        fmt = r[8].split(":")
        val = r[9].split(":")
        if "AF" not in fmt or fmt.index("AF") >= len(val):
            continue
        try:
            afs = [float(x) for x in val[fmt.index("AF")].split(",") if x not in ("", ".")]
            af = max(afs)
        except ValueError:
            continue
        # Heteroplasmic: 5% to 95%. Below 5% NUMT reads and noise dominate.
        if not 0.05 <= af < 0.95:
            continue
        alts = r[4].split(",")
        if any(in_polyc(int(r[1]), r[3], a) for a in alts):
            polyc.append({"variant": f"{r[0]}:{r[1]} {r[3]}>{r[4]}", "af": round(af, 3)})
        else:
            het += 1
            # Each counted site, for the report's chart (bin/report_charts.py).
            sites.append({"pos": int(r[1]), "ref": r[3], "alt": r[4], "af": round(af, 3)})
    ann = mito_annotation(d, s, p)
    listed, plp = [], 0
    for r in ann.pop("rows", []):
        terms = clinvar_terms(r.get("clinvar_significance"))
        is_plp = any(t in CLINVAR_PLP for t in terms)
        plp += is_plp
        try:
            het_row = 0.05 <= float(r.get("af")) < 0.95
        except (TypeError, ValueError):
            het_row = False
        if is_plp or het_row or (terms and not all(t in CLINVAR_BENIGN for t in terms)):
            listed.append({k: r.get(k) for k in MITO_ANNOT_COLS if k not in ("chrom", "filter")}
                          | {"polyc": in_polyc(int(r["pos"]), r["ref"], r["alt"])})
    if "source" in ann:
        ann.update(variants=listed, clinvar_plp=plp)
    return p, {**mito_extended(d, s, p), "pass": passed, "heteroplasmic": het, "heteroplasmic_sites": sites, "heteroplasmy_floor": 0.05,
               "polyc_heteroplasmic": polyc,
               "polyc_tracts": [f"m.{a}-{b}" for a, b in POLYC_TRACTS], "annotation": ann}


def cpsr_cdna_checks(d, s):
    """Step 17's cDNA read check of CPSR's P/LP rows (bin/clip_check.py
    --cpsr), keyed like step 06's; None when step 17 did not write it."""
    p = first_existing(d, [f"cpsr/{s}_cpsr_cdna_check.tsv"])
    return clinvar_cdna_checks(d, s, p) if p else None


def sec_cpsr(d, s):
    html = first_existing(d, [f"cpsr/{s}.cpsr.grch38.html", "cpsr/*.cpsr.*.html"])
    tsv = first_existing(d, [f"cpsr/{s}.cpsr.grch38.classification.tsv.gz", "cpsr/*.cpsr.*.classification.tsv.gz"])
    p = tsv or html
    if not p:
        return None, {}
    out = {"html_report": bool(html)}
    if tsv:
        try:
            with gzip.open(tsv, "rt", errors="replace") as f:
                rows = list(csv.DictReader(f, delimiter="\t"))
                cols = rows[0].keys() if rows else []
        except (OSError, EOFError, gzip.BadGzipFile, zlib.error) as e:
            raise Unreadable(f"{os.path.basename(tsv)}: {e}") from e
        col = next((c for c in ("CLASSIFICATION", "FINAL_CLASSIFICATION", "CPSR_CLASSIFICATION") if c in cols), None)
        if col:
            counts = {}
            for r in rows:
                counts[r.get(col) or "."] = counts.get(r.get(col) or ".", 0) + 1
            out["classification_column"] = col
            out["classification"] = dict(sorted(counts.items(), key=lambda kv: -kv[1]))
            out.update(cpsr_compare(d, s, rows, cols, col))
        else:
            out["classification_column"] = None
    return p, out


def cpsr_compare(d, s, rows, cols, col):
    """CPSR's own class beside the final class and ClinVar's.

    The final class (CLASSIFICATION) is ClinVar's whenever ClinVar has a
    trusted record, so it hides CPSR's own call (CPSR_CLASSIFICATION). This
    lists the P/LP rows (final or own class), counts the final P/LP rows that
    CPSR's own rules class benign or likely benign, and lists every row where
    ClinVar and CPSR's own class disagree on P/LP. ClinVar's class is
    CLINVAR_CLASSIFICATION, or the final class where ASSERTION_AUTHORITY is
    ClinVar. A table with one class column has no disagreement to show.
    Low-confidence: DP_CONTROL below step 6's depth floor (CPSR writes no
    GQ). The cDNA read check of step 17 goes on each P/LP row."""
    own = "CPSR_CLASSIFICATION" if "CPSR_CLASSIFICATION" in cols and col != "CPSR_CLASSIFICATION" else None
    hits = first_existing(d, [f"clinvar/{s}_clinvar_hits.vcf"])
    floor = call_floor(hits) if hits else dict(CALL_FLOOR)
    cdna = cpsr_cdna_checks(d, s)
    norm = clip_check.norm_class
    plp, disagree, own_benign = [], [], 0
    for r in rows:
        final, mine = norm(r.get(col)), norm(r.get(own)) if own else ""
        authority = (r.get("ASSERTION_AUTHORITY") or "").strip()
        clinvar_raw = ((r.get("CLINVAR_CLASSIFICATION") or "").strip() if "CLINVAR_CLASSIFICATION" in cols
                       else (r.get(col) or "").strip() if authority.lower() == "clinvar" else "")
        clinvar = norm(clinvar_raw)
        if final in clip_check.PLP and mine in ("benign", "likely benign"):
            own_benign += 1
        is_plp = final in clip_check.PLP or mine in clip_check.PLP
        differ = bool(own and clinvar and mine) and (clinvar in clip_check.PLP) != (mine in clip_check.PLP)
        if not (is_plp or differ):
            continue
        v = clip_check.cpsr_variant(r)
        dp = (r.get("DP_CONTROL") or "").strip()
        dp = int(dp) if dp.isdigit() else None
        x = {"variant": f"{v[0]}:{v[1]} {v[2]}>{v[3]}" if v else (r.get("GENOMIC_CHANGE") or "."),
             "gene": r.get("SYMBOL") or ".", "final": r.get(col) or ".",
             "cpsr": (r.get(own) or ".") if own else None, "clinvar": clinvar_raw if clinvar else None,
             "authority": authority or None, "acmg_codes": (r.get("ACMG_CODE") or "").strip() or None,
             "dp": dp, "low_confidence": dp is not None and dp < floor["DP"]}
        if cdna is not None and is_plp:
            x["cdna"] = (cdna.get((v[0], v[1], v[2], v[3])) if v else None) or {
                "check": "not_checked", "note": "not in the read-check table", "alt_reads": None,
                "alt_clipped": None, "ref_reads": None, "ref_clipped": None}
        if is_plp:
            plp.append(x)
        if differ:
            disagree.append(x)
    out = {"plp": plp, "disagreements": disagree, "final_plp_own_benign": own_benign,
           "own_class_column": own, "low_confidence": sum(1 for x in plp if x["low_confidence"]),
           "call_floor_dp": floor["DP"]}
    if cdna is not None:
        out["cdna_flagged"] = sum(1 for x in plp if x["cdna"]["check"] == "flagged")
    return out


def csq_rank(fmt, csq):
    """The rule of bin/csq_rank.awk on one CSQ value: (class, entries, freq).
    class is MANE (entries of the MANE Select transcript), else CANONICAL,
    else ANY (every entry); freq is rare, common, unknown (no value) or none
    (the CSQ has no frequency field): MAX_AF, else gnomADe_AF and gnomADg_AF,
    the highest over the entries."""
    ix = {k: i for i, k in enumerate(fmt)}
    entries = [v for v in (t.split("|") for t in csq.split(",") if t) if len(v) == len(fmt)]

    def get(v, k):
        return v[ix[k]] if k in ix else ""
    mane = [v for v in entries if get(v, "MANE_SELECT")]
    canon = [v for v in entries if get(v, "CANONICAL") == "YES"]
    cls, chosen = ("MANE", mane) if mane else ("CANONICAL", canon) if canon else ("ANY", entries)
    fields = ["MAX_AF"] if "MAX_AF" in ix else [f for f in ("gnomADe_AF", "gnomADg_AF") if f in ix]
    afs = []
    for v in entries:
        for f in fields:
            for a in get(v, f).split("&"):
                try:
                    afs.append(float(a))
                except ValueError:
                    pass
    freq = "none" if not fields else "unknown" if not afs else "rare" if max(afs) < COMMON_AF else "common"
    return cls, [{"impact": get(v, "IMPACT"), "gene": get(v, "SYMBOL"),
                  "consequence": get(v, "Consequence").replace("&", ",")} for v in chosen], freq


def call_check(fmt, sample, floor):
    """bin/csq_rank.awk's genotype check of one sample: (weak, reasons, gq,
    dp, ad). Weak: GQ or DP below the floor, a heterozygote's alt read
    fraction outside HET_AB, a homozygous alt call's below HOM_AB."""
    f = dict(zip(fmt.split(":"), sample.split(":")))
    gq, dp = (int(f[k]) if f.get(k, "").isdigit() else None for k in ("GQ", "DP"))
    ad = f.get("AD", "")
    why = []
    if gq is not None and gq < floor["GQ"]:
        why.append(f"GQ {gq} < {floor['GQ']}")
    if dp is not None and dp < floor["DP"]:
        why.append(f"DP {dp} < {floor['DP']}")
    al = f.get("GT", "").replace("|", "/").split("/")
    depths = ad.split(",")
    if len(al) == 2 and all(a.isdigit() for a in al) and len(depths) > 1 and all(x.isdigit() for x in depths):
        tot = sum(int(x) for x in depths)
        if tot:
            ab = sum(int(x) for x in depths[1:]) / tot
            lo, hi = floor["HET_AB"]
            if (al[0] == "0") != (al[1] == "0") and not lo <= ab <= hi:
                why.append(f"alt reads {ab:.2f} of the depth on a heterozygote")
            elif al[0] == al[1] != "0" and ab < floor["HOM_AB"]:
                why.append(f"alt reads {ab:.2f} of the depth on a homozygous call")
    return bool(why), why, gq, dp, ad or None


def clinical_call_floor(path):
    """The floors bin/csq_rank.awk wrote into the clinical VCF, the defaults without the line."""
    floor = dict(CALL_FLOOR, **AB_FLOOR)
    try:
        with open_text(path) as f:
            for line in f:
                if not line.startswith("##"):
                    break
                if line.startswith("##PGP_CALL_FLOOR="):
                    for k, v in re.findall(r"\b(GQ|DP)=(\d+)", line):
                        floor[k] = int(v)
                    m = re.search(r"HET_AB=([0-9.]+)-([0-9.]+)", line)
                    if m:
                        floor["HET_AB"] = (float(m.group(1)), float(m.group(2)))
                    m = re.search(r"HOM_AB=([0-9.]+)", line)
                    if m:
                        floor["HOM_AB"] = float(m.group(1))
    except (OSError, EOFError, gzip.BadGzipFile, zlib.error) as e:
        raise Unreadable(f"{os.path.basename(path)}: {e}") from e
    return floor


def csq_high_impact_genes(path):
    """{(chrom, pos, ref, alt): record} of the records whose ranked CSQ entries
    (bin/csq_rank.awk's rule: the MANE Select transcript, else the canonical
    one, else any) hold a HIGH-impact consequence. record: genotype, genes
    ({gene: consequence}), transcript (MANE, CANONICAL or ANY), frequency
    (rare, common, unknown or none), weak (the call check failed), why, gq,
    dp and ad. A HIGH consequence on a transcript other than the MANE Select
    one does not count when the record has a MANE Select transcript."""
    floor = clinical_call_floor(path)
    fmt, out = None, {}
    try:
        with open_text(path) as f:
            for line in f:
                if line.startswith("##INFO=<ID=CSQ"):
                    m = re.search(r"Format: ([^\"]+)", line)
                    fmt = m.group(1).split("|") if m else None
                    continue
                if line.startswith("#") or not fmt or "IMPACT" not in fmt or "SYMBOL" not in fmt:
                    continue
                r = line.rstrip("\n").split("\t")
                if len(r) < 8:
                    continue
                cls, entries, freq = csq_rank(fmt, info_map(r[7]).get("CSQ", ""))
                genes = {}
                for e in entries:
                    if e["impact"] == "HIGH":
                        genes.setdefault(e["gene"], e["consequence"])
                if genes:
                    gt = r[9].split(":")[0] if len(r) > 9 else ""
                    weak, why, gq, dp, ad = call_check(r[8], r[9], floor) if len(r) > 9 else (False, [], None, None, None)
                    out[(r[0], int(r[1]), r[3], r[4])] = {
                        "genotype": gt, "ps": gt_field(r, "PS"), "genes": genes, "transcript": cls, "frequency": freq,
                        "weak": weak, "why": why, "gq": gq, "dp": dp, "ad": ad}
    except (OSError, EOFError, gzip.BadGzipFile, zlib.error) as e:
        raise Unreadable(f"{os.path.basename(path)}: {e}") from e
    return out


# --- the ACMG SF rules and the carrier list (assets/) --------------------------

_ASSETS = {}


def asset_rows(name):
    """The rows of assets/<name> (comment lines left out), read once."""
    p = os.path.join(ASSETS_DIR, name)
    if p in _ASSETS:
        return _ASSETS[p]
    if not os.path.isfile(p):
        raise Unreadable(f"{name} is not in {ASSETS_DIR}: the report needs the repository's assets/ folder "
                         "(render_report.py --assets)")
    with open(p, encoding="utf-8") as f:
        rows = list(csv.DictReader((x for x in f if not x.startswith("#")), delimiter="\t", quoting=csv.QUOTE_NONE))
    _ASSETS[p] = rows
    return rows


def allele_key(locus):
    """'chr6:26092913:G:A' -> ('chr6', 26092913, 'G', 'A')."""
    c, pos, ref, alt = locus.split(":")
    return (c, int(pos), ref, alt)


def parse_alleles(text, kinds=False):
    """The variants or severity cell of assets/acmg_sf.tsv and
    assets/carrier_genes.tsv: '.' or entries separated by ';', each
    '[kind ]chrom:pos:ref:alt name'. {key: {'kind', 'name'}}."""
    out = {}
    for e in (text or ".").split(";"):
        e = e.strip()
        if e in ("", "."):
            continue
        parts = e.split(" ", 2 if kinds else 1)
        if kinds and len(parts) < 2:
            raise ValueError(f"severity entry needs a kind and a chrom:pos:ref:alt: {e}")
        kind = parts.pop(0) if kinds else None
        if kinds and kind not in SEVERITY_KINDS:
            raise ValueError(f"severity '{kind}' is not one of {', '.join(SEVERITY_KINDS)}: {e}")
        out[allele_key(parts[0])] = {"kind": kind, "name": parts[1] if len(parts) > 1 else parts[0]}
    return out


def acmg_sf_rules():
    """{gene: rule} of assets/acmg_sf.tsv."""
    rules = {}
    for n, r in enumerate(asset_rows("acmg_sf.tsv"), 1):
        try:
            g = r["gene"]
            if r["report_if"] not in SF_REPORT_IF or r["sv_cn"] not in ("yes", "no"):
                raise ValueError(f"report_if must be one of {', '.join(SF_REPORT_IF)}, sv_cn yes or no")
            if r["sv_cn"] == "yes" and r["report_if"] not in ("any", "het_female_too"):
                raise ValueError("sv_cn yes needs a rule that one variant meets (any, het_female_too)")
            regions = []
            for e in (r["sv_regions"] or ".").split(";"):
                e = e.strip()
                if e not in ("", "."):
                    loc, _, name = e.partition(" ")
                    c, _, span = loc.partition(":")
                    a, _, b = span.partition("-")
                    regions.append((c, int(a), int(b), name or loc))
            variants = parse_alleles(r["variants"])
            if r["report_if"] == "hom_variant" and not variants:
                raise ValueError("report_if hom_variant needs its variants")
            if g in rules:
                raise ValueError("the gene is on two rows")
            rules[g] = {"inheritance": r["inheritance"], "report_if": r["report_if"], "variants": variants,
                        "sv_cn": r["sv_cn"] == "yes", "sv_regions": regions,
                        "severity": parse_alleles(r["severity"], kinds=True),
                        "table_text": r["table_text"], "note": None if r["note"] in ("", ".") else r["note"]}
        except (KeyError, ValueError, TypeError) as e:
            raise Unreadable(f"acmg_sf.tsv row {n}: {e}") from e
    if not rules:
        raise Unreadable("acmg_sf.tsv has no gene")
    return rules


def carrier_genes():
    """{gene: row} of assets/carrier_genes.tsv, the severity cell parsed."""
    out = {}
    for n, r in enumerate(asset_rows("carrier_genes.tsv"), 1):
        try:
            out[r["gene"]] = {"inheritance": r["inheritance"], "lists": r["lists"], "condition": r["condition"],
                              "severity": parse_alleles(r.get("severity"), kinds=True),
                              "span": (r["chrom"], int(r["start"]), int(r["end"]))}
        except (KeyError, ValueError, TypeError) as e:
            raise Unreadable(f"carrier_genes.tsv row {n}: {e}") from e
    return out


def allele_copies(gt):
    """(alt copies, ploidy) of a GT: '0/1' (1, 2), '1/1' (2, 2), '1' (1, 1),
    '1/.' (1, 2): a missing allele is a copy not called, not a missing copy,
    so the half of a split 1/2 call is not read as hemizygous."""
    a = [x for x in re.split(r"[/|]", gt or "") if x != ""]
    return sum(1 for x in a if x not in ("0", ".")), len(a)


def zygosity_word(copies, ploidy):
    if not copies:
        return "unknown"
    if ploidy == 1:
        return "hemizygous"
    return "hom" if copies >= 2 else "het"


# --- two hits in one gene across callers: a deletion plus a small variant --------
#
# Small variants alone miss a gene hit twice when one hit is a deletion: an
# exon deletion of one copy of MUTYH plus a pathogenic variant on the other, an
# HBA deletion plus an HBA point variant, one SMN1 copy plus an SMN1 variant, a
# 1q21.1 deletion plus the RBM8A low-expression allele (TAR syndrome: Albers et
# al., Nat Genet 2012;44(4):435-439, PMID 22366785). The deletions come from the
# engines that call them: AnnotSV's class 4 and 5 deletions (step 05), step 35's
# interval table (bin/cn_loci.tsv), its HBA count and Parascopy's SMN1 count.
# Two hits are in trans (biallelic) only when the reads show it: the variant sits
# inside the deletion and is called homozygous (only one copy is left, and it
# carries the variant), or both calls carry a read-backed phase (GT with "|" and
# the same PS) on different copies. Otherwise they are possibly in trans; on the
# same copy (cis) they are no second hit.
TWO_HIT_ORDER = ("biallelic", "possibly biallelic")


def gene_deletions(d, s, spans):
    """The genes a called deletion reaches: ({gene: [deletion]}, engines read,
    engines not read). spans: {gene: (chrom, start, end)}, 1-based inclusive
    (assets/carrier_genes.tsv). A deletion: chrom, start, end (1-based,
    inclusive), text, inside (a variant inside its span is on the one copy left;
    False for the HBA count, which cannot say which gene of the pair is gone),
    gt and ps (its phase, AnnotSV only)."""
    out, read, unread = {}, [], []

    def add(g, x):
        out.setdefault(g, []).append(x)

    def over(chrom, a, b):
        return [g for g, (c, gs, ge) in spans.items() if c == chrom and gs <= b and ge >= a]
    try:
        p, data = sec_annotsv(d, s)
    except (Unreadable, OSError, ValueError, KeyError) as e:
        p, data = None, {}
        unread.append(f"AnnotSV (step 05): {e}")
    if p:
        read.append("AnnotSV class 4 and 5 deletions (step 05)")
        for x in data["rows"]:
            chrom, _, span = x["sv"].partition(":")
            a, _, b = span.partition("-")
            if x["type"] != "DEL" or x["common"] or not (a.isdigit() and b.isdigit()):
                continue
            # A call genotyped 0/0 (or 0|0) is no deletion in this sample: it pairs with nothing.
            if x.get("gt") and set(x["gt"].replace("|", "/").split("/")) == {"0"}:
                continue
            chrom = chrom if chrom.startswith("chr") else "chr" + chrom
            # SV_start is the base before the deleted bases (annotsv_postfilter.sv_span):
            # they are SV_start+1 to SV_end, 1-based. The text keeps AnnotSV's coordinates.
            first = int(a) + 1
            hit = set(over(chrom, first, int(b)))
            # A listed gene AnnotSV names counts only when its span reaches the deleted bases: one
            # that ends at the padding base loses none. A gene with no span here is taken as named.
            named = {g for g in x["genes"] if g not in spans}
            for g in sorted(hit | named):
                add(g, {"chrom": chrom, "start": first, "end": int(b), "inside": True,
                        "text": f"DEL {chrom}:{a}-{b} (AnnotSV class {x['acmg_class']}, {x['caller']})",
                        "gt": x.get("gt"), "ps": x.get("ps")})
    elif not unread:
        unread.append("AnnotSV (step 05): no table")
    try:
        p, data = sec_parascopy(d, s)
        why = "none"
    except (Unreadable, OSError, ValueError, KeyError) as e:
        p, data, why = None, {}, str(e)
    if not p:
        unread.append(f"step 35's copy-number tables: {why}")
        return out, read, unread
    iv = data.get("cn_intervals")
    if iv is not None:
        read.append("step 35's interval table (read depth)")
        for r in iv["changes"]:
            if r["change"] != "deletion" or not (r["start"].isdigit() and r["end"].isdigit()):
                continue
            a, b = int(r["start"]) + 1, int(r["end"])  # the table is BED: 0-based start
            for g in over(r["chrom"], a, b):
                add(g, {"chrom": r["chrom"], "start": a, "end": b, "inside": True, "gt": None, "ps": None,
                        "text": f"{r['locus']} deletion, {r['copies']} of {r['expected']} copies "
                                "(read depth, step 35)"})
    else:
        unread.append("step 35's interval table: none")
    hba = data.get("hba") or {}
    n = hba.get("alpha_globin_genes", "")
    if hba.get("status") == "called" and n.isdigit():
        read.append("step 35's HBA1/HBA2 count")
        # 0 genes: no alpha-globin gene is left to carry a variant, so a hit there pairs with nothing.
        if 1 <= int(n) < 4:
            for g in ("HBA1", "HBA2"):
                if g in spans:
                    c, a, b = spans[g]
                    add(g, {"chrom": c, "start": a, "end": b, "inside": False, "gt": None, "ps": None,
                            "text": f"{n} of 4 alpha-globin genes (read depth, step 35)"})
    else:
        unread.append(f"step 35's HBA1/HBA2 count: {hba.get('message') or 'not called'}")
    smn1 = (data.get("smn1") or {}).get("cn")
    if smn1 is not None:
        read.append("Parascopy's SMN1 count (step 35)")
        if smn1 == 1 and "SMN1" in spans:
            c, a, b = spans["SMN1"]
            add("SMN1", {"chrom": c, "start": a, "end": b, "inside": True, "gt": None, "ps": None,
                         "text": "one SMN1 copy (Parascopy, step 35)"})
    else:
        unread.append("Parascopy's SMN1 count: not resolved")
    return out, read, unread


def phase_haps(gt):
    """The copies (haplotype indexes) a phased GT puts the alt allele on, or
    None when the GT is not phased or a copy is not called: '0|1' -> {1},
    '1|0' -> {0}, '0/1' -> None, '.|1' -> None (the missing copy could carry
    the other allele, so no cis or trans can be read from it)."""
    if not gt or "|" not in gt:
        return None
    alleles = gt.split("|")
    if "." in alleles:
        return None
    return frozenset(i for i, a in enumerate(alleles) if a != "0")


def pair_phase(dl, v):
    """(verdict, why) of a deletion and a small variant in one gene: verdict
    biallelic, possibly biallelic, or cis (no second hit)."""
    inside = dl["inside"] and dl["chrom"] == v["chrom"] and dl["start"] <= v["pos"] <= dl["end"]
    # Read-backed phase first: on the same copy (cis) is no second hit, wherever the variant sits.
    hd, hv = phase_haps(dl.get("gt")), phase_haps(v.get("gt"))
    phased = None
    if hd and hv and dl.get("ps") and dl.get("ps") == v.get("ps"):
        if hd == hv:
            return "cis", f"in cis: read-backed phase puts both on the same copy (phase set {v['ps']})"
        phased = (f"in trans: read-backed phase puts them on different copies (phase set {v['ps']})" if not hd & hv
                  else f"read-backed phase puts the variant on both copies (phase set {v['ps']}); check the reads")
    if inside and v["copies"] >= 2:
        return "biallelic", ("the variant is called homozygous inside the deletion: it is hemizygous, on the one "
                             "copy the deletion leaves")
    if inside and v["ploidy"] == 2:
        return "possibly biallelic", ("a heterozygous call inside a called deletion: one of the two calls may be "
                                      "wrong; check the reads")
    if phased:
        return ("biallelic" if phased.startswith("in trans") else "possibly biallelic"), phased
    return "possibly biallelic", "phase unknown, possibly in trans (no read-backed phase links the two calls)"


def two_hits(variants, dels):
    """The genes hit twice: ({gene: entry}, {(gene, key): deletion text} of
    the variants called homozygous inside a deletion, shown as hemizygous).
    variants: [{gene, key, chrom, pos, text, genotype, copies, ploidy, gt, ps,
    allele}] (allele: the severity name of a hypomorph, or None). An entry is
    the best pair of the gene: verdict, why, deletion, variant, genotype, and
    the allele name; a gene whose every pair is in cis has none. pairs:
    {(gene, key): "<verdict> with <deletion>: <why>"}, the best pair of each
    variant that has one."""
    entries, hemi, pairs = {}, {}, {}
    for v in variants:
        for dl in dels.get(v["gene"], []):
            verdict, why = pair_phase(dl, v)
            if verdict == "biallelic" and why.startswith("the variant is called homozygous"):
                hemi[(v["gene"], v["key"])] = dl["text"]
            if verdict == "cis":
                continue
            e = {"gene": v["gene"], "verdict": verdict, "why": why, "deletion": dl["text"], "variant": v["text"],
                 "genotype": "hemizygous" if (v["gene"], v["key"]) in hemi else v["genotype"],
                 "allele": v.get("allele")}
            old = entries.get(v["gene"])
            if old is None or TWO_HIT_ORDER.index(verdict) < TWO_HIT_ORDER.index(old["verdict"]):
                entries[v["gene"]] = e
            k = (v["gene"], v["key"])
            if k not in pairs or TWO_HIT_ORDER.index(verdict) < TWO_HIT_ORDER.index(pairs[k][0]):
                pairs[k] = (verdict, f"{verdict} with {dl['text']}: {why}")
    return entries, hemi, {k: t for k, (_, t) in pairs.items()}


def sf_verdicts(gene, rule, alleles):
    """Which of a gene's variants the ACMG rule reports. alleles: {key:
    (copies, ploidy)} of the gene's distinct small variants, structural ones
    keyed by their span. {key: None (reported) or the reason it is not}, and
    the partial-deficiency note of a reported gene with a hypomorph allele."""
    sev = rule["severity"]
    out, kept = {}, {}
    for k, cp in alleles.items():
        if sev.get(k, {}).get("kind") == "pseudodeficiency":
            out[k] = f"{sev[k]['name']}: a pseudodeficiency allele, not counted"
        else:
            kept[k] = cp
    n = sum(c for c, _ in kept.values())
    ri = rule["report_if"]
    if ri in ("any", "het_female_too"):
        ok = {k: True for k in kept}
    elif ri == "biallelic":
        ok = {k: n >= 2 for k in kept}
    elif ri == "hemizygous":
        hemi = any(p == 1 for _, p in kept.values())
        ok = {k: hemi or n >= 2 for k in kept}
    else:  # hom_variant
        ok = {k: k in rule["variants"] and c >= 2 and p == 2 for k, (c, p) in kept.items()}
    for k, good in ok.items():
        if good:
            out[k] = None
        elif ri == "hom_variant":
            out[k] = f"{gene} is reported only for: {rule['table_text']}"
        elif ri == "hemizygous":
            out[k] = f"one heterozygous variant; {gene} is reported for {rule['table_text']}"
        else:
            out[k] = f"one heterozygous variant; {gene} is reported with {rule['table_text']}"
    partial = None
    hypo = [sev[k]["name"] for k, why in out.items() if why is None and sev.get(k, {}).get("kind") == "hypomorph"]
    if hypo:
        partial = (f"partial, not profound: {'; '.join(hypo)} is a hypomorph allele, so the deficiency it causes "
                   "is partial")
    return out, partial


def sf_structural(d, s, rules):
    """Step 05's class 4 and 5 SVs (AnnotSV) in genes whose rule allows them
    (sv_cn yes), and the class 4 or 5 deletions over a gene's sv_regions (the
    3' end of EPCAM for MSH2). Common SVs are left out. (rows, checked, note)."""
    try:
        p, data = sec_annotsv(d, s)
    except (Unreadable, OSError, ValueError, KeyError) as e:
        return [], False, f"the AnnotSV table could not be read: {e}"
    if not p:
        return [], False, None
    out = []
    for x in data["rows"]:
        if x["common"] or str(x["acmg_class"]) not in ("4", "5"):
            continue
        chrom, _, span = x["sv"].partition(":")
        a, _, b = span.partition("-")
        chrom = chrom if chrom.startswith("chr") else "chr" + chrom
        hits = {}
        for g in x["genes"]:
            if g in rules and rules[g]["sv_cn"]:
                hits[g] = None
        for g, rule in rules.items():
            for c, rs, re_, name in rule["sv_regions"]:
                if (x["type"] == "DEL" and c == chrom and a.isdigit() and b.isdigit()
                        and int(a) <= re_ and int(b) >= rs and g not in hits):
                    hits[g] = name
        for g, via in sorted(hits.items()):
            out.append({"gene": g, "variant": f"{x['type']} {chrom}:{a}-{b}", "acmg_class": str(x["acmg_class"]),
                        "caller": x["caller"], "support": x["support"], "via": via,
                        "key": ("sv", chrom, a, b, x["type"])})
    return out, True, None


def acmg_sf_tier(d, s, clinical_vcf):
    """The variants that are ACMG secondary findings by the rule of their gene
    (assets/acmg_sf.tsv): the ClinVar P/LP hits of step 6, the clinical
    filter's rare HIGH-impact records, and step 05's class 4 and 5 SVs in the
    genes whose rule allows them. A variant in an SF gene that its rule does
    not report (one MUTYH variant, a heterozygous HFE p.C282Y) is listed
    under not_reported with the reason."""
    rules = acmg_sf_rules()
    clinvar, high = [], []
    hits = first_existing(d, [f"clinvar/{s}_clinvar_hits.vcf"])
    if hits:
        floor = call_floor(hits)
        for r in vcf_records(hits):
            if len(r) < 8:
                continue
            row = clinvar_row(r, floor)
            copies, ploidy = allele_copies(r[9].split(":")[0] if len(r) > 9 else "")
            for g in row["gene"].split(","):
                if g in rules:
                    clinvar.append({"gene": g, "variant": f"{row['chrom']}:{row['pos']} {row['ref']}>{row['alt']}",
                                    "genotype": row["genotype"], "significance": row["significance"],
                                    "stars": row["stars"], "low_confidence": row["low_confidence"],
                                    "key": (row["chrom"], row["pos"], row["ref"], row["alt"]),
                                    "copies": copies, "ploidy": ploidy,
                                    "gt": r[9].split(":")[0] if len(r) > 9 else "", "ps": gt_field(r, "PS")})
    if clinical_vcf:
        # A common record is in the clinical VCF through the ClinVar tier,
        # not as a rare HIGH impact; one of unknown frequency is listed and
        # marked, never counted as rare; a weak call is listed and marked.
        for (c, pos, ref, alt), x in sorted(csq_high_impact_genes(clinical_vcf).items()):
            if x["frequency"] == "common":
                continue
            for g, cons in sorted(x["genes"].items()):
                if g in rules:
                    gt = x["genotype"]
                    copies, ploidy = allele_copies(gt)
                    zyg = {"0/1": "het", "1/0": "het", "0|1": "het", "1|0": "het", "1/1": "hom", "1|1": "hom"}.get(gt, gt)
                    high.append({"gene": g, "variant": f"{c}:{pos} {ref}>{alt}", "genotype": zyg, "consequence": cons,
                                 "transcript": x["transcript"], "frequency": x["frequency"],
                                 "frequency_unknown": x["frequency"] == "unknown",
                                 "weak_call": x["weak"], "call_note": "; ".join(x["why"]) or None,
                                 "gq": x["gq"], "dp": x["dp"], "ad": x["ad"],
                                 "key": (c, pos, ref, alt), "copies": copies, "ploidy": ploidy,
                                 "gt": gt, "ps": x.get("ps")})
    svs, sv_checked, sv_note = sf_structural(d, s, rules)
    # One verdict per gene, over its distinct variants: a ClinVar hit that is
    # also a HIGH-impact record is one variant.
    alleles = {}
    for x in clinvar + high:
        alleles.setdefault(x["gene"], {}).setdefault(x["key"], (max(x["copies"], 1), x["ploidy"] or 2))
    for x in svs:
        alleles.setdefault(x["gene"], {}).setdefault(x["key"], (1, 2))
    # A called deletion of a gene whose rule needs two variants counts as one
    # of them (two_hits), unless read-backed phase puts it on the variant's copy.
    two_hit_note = None
    try:
        spans = {g: x["span"] for g, x in carrier_genes().items()}
    except Unreadable as e:
        spans = {}
        two_hit_note = (f"A deletion as the second hit: only AnnotSV's gene names were read; step 35's deletions "
                        f"need the gene spans of {e}")
    dels, _, _ = gene_deletions(d, s, spans)
    seen, variants = set(), []
    for x in clinvar + high:
        pseudo = rules[x["gene"]]["severity"].get(x["key"], {}).get("kind") == "pseudodeficiency"
        if (x["gene"], x["key"]) not in seen and not pseudo:
            seen.add((x["gene"], x["key"]))
            variants.append({"gene": x["gene"], "key": x["key"], "chrom": x["key"][0], "pos": x["key"][1],
                             "text": x["variant"], "genotype": x["genotype"], "copies": x["copies"],
                             "ploidy": x["ploidy"] or 2, "gt": x["gt"], "ps": x["ps"], "allele": None})
    entries, hemi, pairs = two_hits(variants, dels)
    two = {}
    for g in entries:
        if rules[g]["report_if"] in ("biallelic", "hemizygous"):
            alleles[g][("deletion", g)] = (1, 2)
            two.update({k: t for k, t in pairs.items() if k[0] == g})
    verdict = {g: sf_verdicts(g, rules[g], a) for g, a in alleles.items()}
    reported = {"clinvar": [], "high": [], "sv": []}
    not_reported = []
    for kind, rows, source in (("clinvar", clinvar, "ClinVar"), ("high", high, "HIGH impact"), ("sv", svs, "AnnotSV")):
        for x in rows:
            why, partial = verdict[x["gene"]][0][x["key"]], verdict[x["gene"]][1]
            rule = rules[x["gene"]]
            key = x["key"]
            if (x["gene"], key) in hemi:
                x["genotype"] = "hemizygous"
            for k in ("key", "copies", "ploidy", "gt", "ps"):
                x.pop(k, None)
            if why is None:
                x.update(rule=rule["table_text"], rule_note=rule["note"], severity=partial)
                if kind != "sv" and (x["gene"], key) in two:
                    x["two_hit"] = two[(x["gene"], key)]
                reported[kind].append(x)
            else:
                not_reported.append({"gene": x["gene"], "variant": x["variant"], "genotype": x.get("genotype"),
                                     "source": source, "reason": why})
    high = reported["high"]
    sure = [x for x in high if not (x["weak_call"] or x["frequency_unknown"])]
    return {"version": ACMG_SF_VERSION, "genes_on_list": len(rules), "rules": "assets/acmg_sf.tsv",
            "clinvar_hits": reported["clinvar"], "clinvar_checked": bool(hits),
            "high_impact": high, "high_impact_checked": bool(clinical_vcf),
            "high_impact_rare": len(sure),
            "high_impact_weak": sum(1 for x in high if x["weak_call"]),
            "high_impact_frequency_unknown": sum(1 for x in high if x["frequency_unknown"]),
            "structural": reported["sv"], "structural_checked": sv_checked, "structural_note": sv_note,
            "not_reported": not_reported, "two_hit_note": two_hit_note}


def gt_field(r, key):
    """One FORMAT field of a VCF record's first sample, or None."""
    if len(r) < 10:
        return None
    v = dict(zip(r[8].split(":"), r[9].split(":"))).get(key)
    return None if v in (None, "", ".") else v


def carrier_list(p, d=None, s=None):
    """The sample's own carrier list: the ClinVar hits of step 6 (p) in the
    genes of assets/carrier_genes.tsv, one person, no partner needed. An
    allele the list marks hypomorph or pseudodeficiency is set aside, never a
    carrier finding on its own. With the sample folder (d, s), the genes of
    autosomal recessive conditions that a called deletion and a hit both reach
    (two_hits), and a hit called homozygous inside a deletion is shown as
    hemizygous."""
    genes = carrier_genes()
    carriers, aside, variants = [], [], []
    for r in vcf_records(p):
        if len(r) < 8:
            continue
        row = clinvar_row(r)
        gt = r[9].split(":")[0] if len(r) > 9 else ""
        copies, ploidy = allele_copies(gt)
        key = (row["chrom"], row["pos"], row["ref"], row["alt"])
        for g in row["gene"].split(","):
            if g not in genes:
                continue
            x = genes[g]
            item = {"gene": g, "variant": f"{row['chrom']}:{row['pos']} {row['ref']}>{row['alt']}",
                    "genotype": row["genotype"] if half_call(gt) else zygosity_word(copies, ploidy),
                    "significance": row["significance"],
                    "stars": row["stars"], "inheritance": x["inheritance"], "condition": x["condition"]}
            sev = x["severity"].get(key)
            if sev:
                item.update(kind=sev["kind"], name=sev["name"])
                aside.append(item)
            else:
                carriers.append(item)
            if not sev or sev["kind"] == "hypomorph":
                variants.append({"gene": g, "key": key, "chrom": row["chrom"], "pos": row["pos"],
                                 "text": item["variant"], "genotype": item["genotype"], "copies": copies,
                                 "ploidy": ploidy, "gt": gt, "ps": gt_field(r, "PS"),
                                 "allele": sev["name"] if sev else None, "item": item})
    out = {"genes_on_list": len(genes), "list": "assets/carrier_genes.tsv",
           "carriers": carriers, "set_aside": aside}
    if d is not None:
        dels, read, unread = gene_deletions(d, s, {g: x["span"] for g, x in genes.items()})
        entries, hemi, _ = two_hits(variants, dels)
        for v in variants:
            if (v["gene"], v["key"]) in hemi:
                v["item"].update(genotype="hemizygous",
                                 hemizygous=f"called homozygous inside {hemi[(v['gene'], v['key'])]}")
        rows = []
        for g, e in sorted(entries.items()):
            if genes[g]["inheritance"] == "AR":
                e.update(condition=genes[g]["condition"], inheritance="AR")
                rows.append(e)
        out["two_hits"] = {"genes": rows, "deletions_read": read, "deletions_not_read": unread}
    return out


def acmg_sf_or_note(d, s, clinical_vcf):
    """acmg_sf_tier, or the reason it could not run (a missing or broken
    assets/acmg_sf.tsv), which both reports print: the clinical counts stay."""
    try:
        return acmg_sf_tier(d, s, clinical_vcf)
    except Unreadable as e:
        return {"version": ACMG_SF_VERSION, "unavailable": str(e), "clinvar_hits": [], "high_impact": [],
                "structural": [], "not_reported": []}


def carrier_list_or_note(p, d=None, s=None):
    """carrier_list, or the reason it could not run, which both reports print:
    the ClinVar hits stay."""
    try:
        return carrier_list(p, d, s)
    except Unreadable as e:
        return {"unavailable": str(e)}


def sec_clinical(d, s):
    p = first_existing(d, [f"clinical/{s}_clinical_summary.tsv"])
    v = first_existing(d, [f"clinical/{s}_clinical.vcf.gz"])
    if not p:
        # The Nextflow report gets the clinical VCF, not the summary table:
        # the count only (one record per variant, as in the table).
        if not v:
            return None, {}
        return v, {"variants": sum(1 for _ in vcf_records(v)), "acmg_sf": acmg_sf_or_note(d, s, v)}
    rows = read_tsv(p)
    by = {}
    for r in rows:
        by[r.get("IMPACT") or "."] = by.get(r.get("IMPACT") or ".", 0) + 1
    genes = sorted({r.get("GENE") for r in rows if r.get("GENE") not in (None, "", ".")})
    out = {"variants": len(rows), "by_impact": by, "genes": len(genes), "acmg_sf": acmg_sf_or_note(d, s, v)}
    # Step 23's tier and call check columns; a table written before them has neither.
    if rows and "FREQ_TIER" in rows[0]:
        out["frequency_unknown"] = sum(1 for r in rows if r.get("FREQ_TIER") == "unknown")
    if rows and "CALL_CHECK" in rows[0]:
        out["weak_calls"] = sum(1 for r in rows if r.get("CALL_CHECK") == "weak")
    return p, out


def sec_slivar(d, s):
    p = first_existing(d, [f"slivar/{s}_slivar_summary.tsv"])
    if not p:
        # The Nextflow report gets the prioritized VCF: the count, plus the
        # phase when COMPHET_PHASE's pairs table sits beside it (HTML_REPORT
        # stages it there).
        v = first_existing(d, [f"slivar/{s}_prioritized.vcf.gz"])
        if not v:
            return None, {}
        out = {"prioritized": sum(1 for _ in vcf_records(v))}
        out.update(comphet_phase(os.path.join(os.path.dirname(v), f"{s}_compound_het_pairs.tsv")))
        return v, out
    rows = read_tsv(p)
    out = {"prioritized": len(rows)}
    # Step 31's frequency tier: prioritized for the compound-het search, not rare.
    if rows and "FREQ_TIER" in rows[0]:
        out["frequency_unknown"] = sum(1 for r in rows if r.get("FREQ_TIER") == "unknown")
    ch = os.path.join(os.path.dirname(p), f"{s}_compound_hets.tsv")
    if os.path.isfile(ch):
        crow = read_tsv(ch)
        out["compound_het_variants"] = len(crow)
        out["compound_het_genes"] = len({r.get("GENE") for r in crow if r.get("GENE")})
    out.update(comphet_phase(os.path.join(os.path.dirname(p), f"{s}_compound_het_pairs.tsv")))
    return p, out


# The trans pairs the report lists by name; the rest are counted.
COMPHET_TRANS_SHOWN = 20


def comphet_phase(path):
    """Step 31's pairs table (bin/comphet_phase.py): the pairs by phase from
    the reads, and the trans ones (both copies of the gene hit) by name. A run
    from before the table has none of these keys."""
    if not os.path.isfile(path):
        return {}
    rows = read_tsv(path)
    phase = {k: sum(1 for r in rows if r.get("PHASE") == k) for k in ("trans", "cis", "unknown")}

    def n(x):
        return int(x) if (x or "").isdigit() else 0
    trans = [{"gene": r.get("GENE", ""),
              "variants": f"{r.get('CHROM')}:{r.get('POS_1')} {r.get('REF_1')}>{r.get('ALT_1')} + "
                          f"{r.get('CHROM')}:{r.get('POS_2')} {r.get('REF_2')}>{r.get('ALT_2')}",
              "reads": f"{r.get('TRANS_READS')} of {n(r.get('TRANS_READS')) + n(r.get('CIS_READS'))}"}
             for r in rows if r.get("PHASE") == "trans"]
    return {"compound_het_pairs": len(rows), "compound_het_phase": phase,
            "compound_het_trans": trans[:COMPHET_TRANS_SHOWN]}


def sec_coverage(d, s):
    p = first_existing(d, [f"mosdepth/{s}.mosdepth.summary.txt", f"coverage/{s}.mosdepth.summary.txt"])
    if not p:
        return None, {}
    mean = None
    for r in read_tsv(p):
        if r.get("chrom") == "total":
            mean = r.get("mean")
    return p, {"mean_depth": float(mean) if mean not in (None, "") else None}


# A gene of assets/callability_genes.tsv counts as read when at least this
# share of its bases has the depth step 16b counts (its thresholds column,
# 10x of reads with MAPQ 20 or more).
CALLABLE_SHARE = 0.9


def callability(catalog, table):
    """The reading of every gene of the catalog rows from step 16b's per-gene
    table (mosdepth --thresholds over the catalog's genes: chrom, start, end,
    gene, bases at the threshold), as the section's data."""
    with open_text(table) as f:
        lines = [x.rstrip("\n").split("\t") for x in f if x.strip()]
    if not lines or not lines[0][0].startswith("#") or len(lines[0]) != 5 \
            or not re.fullmatch(r"[0-9]+X", lines[0][4]):
        raise Unreadable(f"{os.path.basename(table)} is not mosdepth's thresholds table of one depth")
    min_depth = int(lines[0][4][:-1])
    bases = {}
    for c in lines[1:]:
        if len(c) != 5:
            raise Unreadable(f"{os.path.basename(table)}: a row has {len(c)} columns, not 5")
        bases[c[3]] = (int(c[2]) - int(c[1]), int(c[4]))
    if not bases:
        raise Unreadable(f"{os.path.basename(table)} has no gene rows")
    # mosdepth writes no row on a contig the BAM lacks: a chromosome with no row at all is
    # not on the sample's reference; a gene missing from a chromosome that has rows means
    # the table was written from another catalog.
    contigs = {c[0] for c in lines[1:]}
    genes, unowned, owned, not_on_reference = [], [], [], []
    for r in catalog:
        g = {"gene": r["gene"], "hgnc_id": r["hgnc_id"], "reasons": r["reasons"].split(","),
             "owner": None if r["owner"] == "." else r["owner"], "share": None}
        if r["start"] == "." or r["chrom"] not in contigs:
            g["status"] = "not on this reference"
            not_on_reference.append(g["gene"])
        else:
            if r["gene"] not in bases:
                raise Unreadable(f"{os.path.basename(table)} has no row for {r['gene']}: it was written from "
                                 "another assets/callability_genes.tsv; run step 16b again")
            n, ok = bases[r["gene"]]
            g["share"] = round(ok / n, 4) if n else 0.0
            g["status"] = "read" if g["share"] >= CALLABLE_SHARE else "not read"
        genes.append(g)
        if g["status"] != "read":
            (owned if g["owner"] else unowned).append(g)
    return {"min_depth": min_depth, "min_mapq": 20, "min_share": CALLABLE_SHARE, "catalog_genes": len(genes),
            "read": sum(1 for g in genes if g["status"] == "read"),
            "not_read": [{k: g[k] for k in ("gene", "hgnc_id", "reasons", "share", "status")} for g in unowned],
            "not_read_owned": [{k: g[k] for k in ("gene", "owner", "share", "status")} for g in owned],
            "not_on_reference": not_on_reference}


def sec_callability(d, s):
    """Step 16b's gene callability (the MOSDEPTH task's in a Nextflow folder)."""
    p = first_existing(d, [f"mosdepth/{s}_callability.bed.gz", f"coverage/{s}_callability.bed.gz"])
    if not p:
        return None, {}
    return p, callability(asset_rows("callability_genes.tsv"), p)


# How far a chrX or chrY copy number may sit from a whole number and still be
# read as that number. Further off (a mosaic, a mixed sample, a noisy run)
# names no karyotype.
KARYOTYPE_TOL = 0.3


# The sex-chromosome aneuploidies named, by (chrX, chrY) copies. Any other pair of
# whole numbers (chrY 4, chrX 6, a degenerate 100000) names nothing: in a run it is
# far more often depth noise or a broken input than a karyotype.
KARYOTYPES = {(1, 0): "45,X", (2, 1): "47,XXY", (1, 2): "47,XYY", (3, 0): "47,XXX",
              (2, 2): "48,XXYY", (3, 1): "48,XXXY", (4, 0): "48,XXXX",
              (4, 1): "49,XXXXY", (5, 0): "49,XXXXX"}


def karyotype(cn_x, cn_y, tol=KARYOTYPE_TOL):
    """The sex-chromosome karyotype the two copy numbers name, one of
    KARYOTYPES; None for 46,XX and 46,XY, for a missing value, for a value
    between whole numbers, and for any pair of copies not in KARYOTYPES."""
    x, y = num(cn_x), num(cn_y)
    if x is None or y is None or not (-1 < x < 10 and -1 < y < 10):   # also keeps inf from round()
        return None
    rx, ry = int(round(x)), int(round(y))
    if abs(x - rx) > tol or abs(y - ry) > tol:
        return None
    return KARYOTYPES.get((rx, ry))


def sec_sex_check(d, s):
    p = first_existing(d, ["indexcov/indexcov-indexcov.ped", "indexcov/*-indexcov.ped"])
    if not p:
        return None, {}
    with open(p, errors="replace") as f:
        lines = [l.rstrip("\n").split("\t") for l in f if l.strip()]
    head = [h.lstrip("#") for h in lines[0]] if lines else []
    if "sex" not in head or len(lines) < 2:
        raise Unreadable(f"{os.path.basename(p)} has no sex column or no sample row")
    row = dict(zip(head, lines[-1]))
    sex = {"1": "male", "2": "female"}.get(row.get("sex", ""), "unknown")
    return p, {"inferred_sex": sex, "cn_chrX": row.get("CNchrX", ""), "cn_chrY": row.get("CNchrY", ""),
               "karyotype": karyotype(row.get("CNchrX"), row.get("CNchrY"))}


# --- sample identity and contamination (step 33) ----------------------------------

FREEMIX_WARN = 0.03   # VerifyBamID2 FREEMIX above this is reported as possible contamination
SAME_PERSON = 0.9     # somalier relatedness at or above this: the same person (a duplicate or a swap)
SOMALIER_SEX = {"1": "male", "2": "female"}   # its samples.tsv sex column; -9 unknown, -2 X and Y disagree


def read_rows(path):
    """Rows of a tab-separated file whose header may start with '#'."""
    with open_text(path) as f:
        lines = [l.rstrip("\n").split("\t") for l in f if l.strip()]
    if not lines:
        return []
    head = [lines[0][0].lstrip("#")] + lines[0][1:]
    return [dict(zip(head, r)) for r in lines[1:]]


def num(x):
    try:
        v = float(x)
    except (TypeError, ValueError):
        return None
    return v if v == v else None   # NaN is no value


def somalier_x_sex(row):
    """male or female when somalier's own rule calls the sex from chrX, else
    None. The rule of relate --infer (relate.nim, add_parents_and_check_sex):
    more than 10 chrX sites with depth, heterozygous / homozygous-ALT sites
    below 0.05 for male or above 0.4 for female, and fewer than 6% of all
    sites with an allele balance outside 0.1 to 0.9."""
    n, het, hom = num(row.get("X_n")), num(row.get("X_het")), num(row.get("X_hom_alt"))
    # A missing column passes (older somalier, test stubs); a present value
    # that is not a number (nan from 0/0 when no autosomal site has a call)
    # fails, as somalier's own `< 0.06` does.
    raw_mid = row.get("p_middling_ab")
    mid = num(raw_mid)
    if n is None or het is None or hom is None or n <= 10 or (
            raw_mid is not None and (mid is None or mid >= 0.06)):
        return None
    if hom == 0:
        return "female" if het > 0 else None
    ratio = het / hom
    return "male" if ratio < 0.05 else "female" if ratio > 0.4 else None


def sample_qc_table(sample, samples_tsv, selfsm=None, pairs_tsv=None, declared_sex=None,
                    freemix_warn=FREEMIX_WARN, somalier_id=None, marker_check=""):
    """The step 33 verdict for one sample, as ordered (key, value) pairs."""
    sid = somalier_id or sample
    rows = read_rows(samples_tsv)
    row = next((r for r in rows if r.get("sample_id") == sid), None)
    if row is None:
        raise Unreadable(f"{os.path.basename(samples_tsv)} has no row for sample '{sid}' "
                         f"(it has: {', '.join(r.get('sample_id', '?') for r in rows) or 'none'})")
    inferred = SOMALIER_SEX.get(row.get("sex", ""), "unknown")
    declared = (declared_sex or "").lower() or None
    # With the declared sex as its pedigree (--ped), somalier starts its sex
    # column from that sex and overwrites it only when the reads tell, so the
    # pedigree's sex left standing is not a call: unknown, as without a
    # pedigree. -2 from a pedigree female with chrY reads is not one either
    # unless chrX looks female too.
    ped = {"male": "male", "female": "female"}.get(row.get("original_pedigree_sex", "").lower())
    x_says = somalier_x_sex(row) if ped else None
    if ped and inferred == ped and x_says != ped:
        inferred = "unknown"
    x_and_y = row.get("sex", "") == "-2" and (not ped or x_says == "female")
    # -2 first: it is worth saying even when no sex was declared
    if x_and_y:
        sex_check = "not_checked"
        why = ("chrX is heterozygous like a female sample but chrY has reads: "
               "a sex-chromosome aneuploidy or a mixed sample")
    elif not declared:
        sex_check, why = "not_checked", "no declared sex"
    elif inferred == "unknown":
        sex_check = "not_checked"
        why = (f"somalier could not tell the sex from {row.get('X_n', '0')} chrX sites "
               f"(it needs more than 10 with reads, and allele balances that look like one person)")
    elif inferred == declared:
        sex_check, why = "ok", f"declared {declared}, somalier infers {inferred}"
    else:
        sex_check, why = "mismatch", f"declared {declared}, somalier infers {inferred}"
    gt = num(row.get("gt_depth_mean"))
    ydp = num(row.get("Y_depth_mean"))
    xdp = num(row.get("X_depth_mean"))
    sites = sum(int(num(row.get(k)) or 0) for k in ("n_hom_ref", "n_het", "n_hom_alt"))
    out = [("sample", sample), ("somalier_id", sid), ("declared_sex", declared or ""),
           ("inferred_sex", inferred), ("sex_check", sex_check), ("sex_check_reason", why),
           ("sites_genotyped", str(sites)), ("depth_at_sites", row.get("gt_depth_mean", "")),
           ("x_sites", row.get("X_n", "")), ("x_het", row.get("X_het", "")),
           ("x_hom_alt", row.get("X_hom_alt", "")),
           ("x_depth_ratio", f"{xdp / gt:.2f}" if gt and xdp is not None else ""),
           ("y_sites", row.get("Y_n", "")),
           ("y_depth_ratio", f"{ydp / gt:.2f}" if gt and ydp is not None else "")]
    freemix, status, markers = None, "not_run", ""
    if selfsm:
        sm = read_rows(selfsm)
        freemix = num(sm[0].get("FREEMIX")) if sm else None
        if freemix is None:
            raise Unreadable(f"{os.path.basename(selfsm)} has no FREEMIX value")
        status = "warn" if freemix > freemix_warn else "ok"
        # #SNPS of the selfSM is the size of the panel, not the markers with reads
        markers = sm[0].get("#SNPS", "")
    # marker_check: passed, or skipped when fewer than 1,000 panel markers had
    # reads and VerifyBamID2 ran with --DisableSanityCheck.
    out += [("freemix", f"{freemix:.4f}" if freemix is not None else ""),
            ("freemix_warn_above", f"{freemix_warn:g}"), ("contamination", status),
            ("panel_markers", markers), ("verifybamid2_marker_check", marker_check)]
    same = []
    if pairs_tsv:
        for r in read_rows(pairs_tsv):
            a, b = r.get("sample_a"), r.get("sample_b")
            rel = num(r.get("relatedness"))
            if sid in (a, b) and a != b and rel is not None and rel >= SAME_PERSON:
                same.append(f"{b if a == sid else a} ({rel:.2f})")
    out.append(("same_person_as", ", ".join(same)))
    return out


def sec_sample_qc(d, s):
    p = first_existing(d, [f"qc/{s}_sample_qc.tsv"])
    if not p:
        return None, {}
    data = {}
    for r in read_tsv(p):
        if r.get("key"):
            data[r["key"]] = r.get("value", "")
    if "inferred_sex" not in data:
        raise Unreadable(f"{os.path.basename(p)} has no inferred_sex row")
    # somalier's chrX and chrY depth over the autosomal depth: one copy reads 0.5.
    # Named only from more than 10 chrX sites, the floor somalier's own sex
    # call needs (somalier_x_sex); a thinner ratio is noise.
    xr, yr = num(data.get("x_depth_ratio")), num(data.get("y_depth_ratio"))
    xs = num(data.get("x_sites"))
    data["karyotype"] = (karyotype(2 * xr, 2 * yr)
                         if xr is not None and yr is not None and xs is not None and math.isfinite(xs) and xs > 10
                         else None)
    return p, data


# --- mosaic loss of chrY or of one chrX (step 33, bin/mosaic_sex_loss.py) ----------

MOSAIC_CALLS = ("detected", "not_detected", "gain", "not_estimated")
MOSAIC_NUMBERS = ("depth_ratio", "depth_ratio_expected", "loss_depth", "loss_depth_low", "loss_depth_high",
                  "loss_af", "loss_af_low", "loss_af_high", "floor")


def sec_mosaic_sex_loss(d, s):
    """The mosaic_* rows of the step 33 table, without the prefix; none in a
    table written before the estimate existed."""
    p = first_existing(d, [f"qc/{s}_sample_qc.tsv"])
    if not p:
        return None, {}
    data = {r["key"][len("mosaic_"):]: r.get("value") or "" for r in read_tsv(p)
            if (r.get("key") or "").startswith("mosaic_")}
    if not data:
        return None, {}
    if data.get("call") not in MOSAIC_CALLS:
        raise Unreadable(f"{os.path.basename(p)} has no mosaic_call row ({', '.join(MOSAIC_CALLS)})")
    for k in MOSAIC_NUMBERS:
        v = num(data.get(k))
        data[k] = v if v is not None and abs(v) != float("inf") else None   # inf is no value either
    data["x_het_sites"] = int(data["x_het_sites"]) if (data.get("x_het_sites") or "").isdigit() else None
    return p, data


def sec_read_count(d, s):
    """Step 02's read count (and the ALIGN_MINIMAP2 task's): the reads that
    went into the aligner and the primary records of the BAM it wrote."""
    p = first_existing(d, [f"qc/{s}_read_count.tsv"])
    if not p:
        return None, {}
    data = {}
    for r in read_tsv(p):
        if r.get("key"):
            data[r["key"]] = r.get("value", "")
    for k in ("input_reads", "primary_records", "missing_reads"):
        if not re.fullmatch(r"-?[0-9]+", data.get(k, "")):
            raise Unreadable(f"{os.path.basename(p)} has no whole number for {k}")
        data[k] = int(data[k])
    if data.get("read_count_check") not in ("ok", "mismatch"):
        raise Unreadable(f"{os.path.basename(p)} has no read_count_check row (ok or mismatch)")
    return p, data


# Sections that are not a step of their own, or not one of one person (the
# couple screen needs two): no "not run" line without them.
NOT_A_STEP = {"read_count", "couple", "carrier"}

SECTIONS = [
    # key, title, step, reader
    ("coverage", "Coverage (mosdepth)", "16b", sec_coverage),
    ("callability", "Gene callability (mosdepth)", "16b", sec_callability),
    ("sex_check", "Sex check (indexcov)", "16", sec_sex_check),
    ("sample_qc", "Sample identity and contamination (somalier, VerifyBamID2)", "33", sec_sample_qc),
    ("mosaic_sex_loss", "Mosaic loss of chromosome Y or X", "33", sec_mosaic_sex_loss),
    ("read_count", "Reads kept by alignment", "02", sec_read_count),
    ("variants", "Variant calling", "03", sec_variants),
    ("clinvar", "ClinVar screen", "06", sec_clinvar),
    ("couple", "Couple carrier screen", "39", sec_couple),
    ("carrier", "Carrier list", "06", sec_carrier),
    ("pharmcat", "PharmCAT", "07", sec_pharmcat),
    ("cpic", "CPIC drug-gene recommendations", "27", sec_cpic),
    ("pypgx", "pypgx", "32", sec_pypgx),
    ("cyrius", "Cyrius (CYP2D6)", "21", sec_cyrius),
    ("parascopy", "SMN1/SMN2 copy number (Parascopy, short reads)", "35", sec_parascopy),
    ("paralogs", "Paralog genes (Paraphase, PacBio HiFi)", "35b", sec_paralogs),
    ("hla", "HLA typing (T1K)", "08", sec_hla),
    ("prs", "Polygenic risk scores", "25", sec_prs),
    ("ancestry_fine", "Fine-grained ancestry (ADMIXTURE, RFMix)", "26b", sec_ancestry_fine),
    ("manta", "Structural variants (Manta)", "04", sec_manta),
    ("delly", "Structural variants (Delly)", "19", sec_delly),
    ("cnvpytor", "Copy number (CNVpytor)", "18", sec_cnvpytor),
    ("sv_consensus", "SV consensus merge", "22", sec_sv_consensus),
    ("annotsv", "Structural variants, class 4 and 5 (AnnotSV)", "05", sec_annotsv),
    ("expansions", "Repeat expansions (ExpansionHunter)", "09", sec_expansions),
    ("telomere", "Telomere content (TelomereHunter)", "10", sec_telomere),
    ("roh", "Runs of homozygosity", "11", sec_roh),
    ("haplogroup", "Mitochondrial haplogroup", "12", sec_haplogroup),
    ("y_haplogroup", "Y-chromosome haplogroup (Yleaf)", "37", sec_y_haplogroup),
    ("y_tree", "Y placement on the current public tree", "37", sec_y_tree),
    ("y_str", "Y-STR markers (HipSTR)", "45", sec_y_str),
    ("consumer_export", "Raw-data export (23andMe v5 layout)", "40", sec_consumer_export),
    ("bloodtype", "Blood groups (predicted)", "41", sec_bloodtype),
    ("lpa", "Lp(a): LPA KIV-2 repeat count", "46", sec_lpa),
    ("vntyper", "MUC1 VNTR frameshifts (VNtyper)", "52", sec_vntyper),
    ("locityper", "Locityper second opinion at hard loci", "53", sec_locityper),
    ("mito", "Mitochondrial variants", "20", sec_mito),
    ("traits", "Traits panel", "42", sec_traits),
    ("archaic", "Archaic ancestry (IBDmix)", "43", sec_archaic),
    ("ancient", "Ancient populations (AADR PCA)", "44", sec_ancient),
    ("viral", "Viral DNA in the unmapped reads", "47", sec_viral),
    ("exomiser", "Symptom-driven ranking (Exomiser)", "48", sec_exomiser),
    ("chip", "Clonal haematopoiesis (CHIP)", "49", sec_chip),
    ("mei", "Mobile element insertions (MEGAnE)", "50", sec_mei),
    ("mosaic", "Mosaic chromosomal alterations (MoChA)", "51", sec_mocha),
    ("cpsr", "Cancer predisposition (CPSR)", "17", sec_cpsr),
    ("clinical", "Clinical variant filter", "23", sec_clinical),
    ("slivar", "Variant prioritization (slivar)", "31", sec_slivar),
]


# --- run status and manifest -----------------------------------------------------

def read_run_status(d):
    p = os.path.join(d, "logs", "run_status.tsv")
    if not os.path.isfile(p):
        return None
    st = {"started_epoch": None, "started_utc": None, "declared_sex": None, "steps": {}}
    with open(p, errors="replace") as f:
        for line in f:
            c = line.rstrip("\n").split("\t")
            if len(c) < 3 or c[0].startswith("#"):
                continue
            if c[0] == "meta" and c[1] == "started_epoch":
                try:
                    st["started_epoch"] = float(c[2])
                except ValueError:
                    pass
            elif c[0] == "meta":
                st[c[1]] = c[2]
            elif c[0] == "step":
                st["steps"][c[1]] = c[2]
    if st["started_epoch"] is not None:
        st["started_utc"] = iso(st["started_epoch"])
    return st


def read_manifest(d):
    p = os.path.join(d, "run_manifest.tsv")
    if not os.path.isfile(p):
        return None
    rows = []
    with open(p, errors="replace") as f:
        for line in f:
            c = line.rstrip("\n").split("\t")
            if len(c) >= 3 and not c[0].startswith("#"):
                rows.append([c[0], c[1], "\t".join(c[2:])])
    return rows


def manifest_value(rows, section, key):
    for r in rows or []:
        if r[0] == section and r[1] == key:
            return r[2]
    return None


def read_consensus(d, s, sections):
    """Step 36's CYP2D6 row (pgx_consensus/<id>_pgx_consensus.tsv), or None
    when the table is missing or older than a caller's result it was made from."""
    p = first_existing(d, [f"pgx_consensus/{s}_pgx_consensus.tsv"])
    if not p:
        return None
    for k in ("pypgx", "cyrius"):
        src = sections[k]["source"]
        if src and os.path.getmtime(os.path.join(d, src)) > os.path.getmtime(p):
            return None
    row = next((r for r in read_tsv(p) if (r.get("Gene") or "").strip() == "CYP2D6"), None)
    if row is None:
        return None
    return {"result": row.get("Result") or "", "passed_to_pharmcat": row.get("Outside_call") == "yes",
            "reason": row.get("Reason") or "", "source": os.path.relpath(p, d)}


def cyp2d6_block(sample, sample_dir, sections):
    """CYP2D6 from each caller, side by side, judged by step 36's rule
    (bin/pgx_outside_calls.py): a no-call in any case (Indeterminate, None,
    '.'...) is not a call, and Cyrius counts only with Filter PASS. PharmCAT
    calls no CYP2D6 from a VCF: its column is what step 36 passed to it, so
    agreement is between pypgx and Cyrius. With step 36's table its verdict
    is kept as well."""
    pc = next((g for g in sections["cpic"]["data"].get("genes", []) if g["gene"] == "CYP2D6"), None)
    pg = sections["pypgx"]["data"].get("cyp2d6") if sections["pypgx"]["state"] == "ok" else None
    cy = sections["cyrius"]["data"] if sections["cyrius"]["state"] == "ok" else None
    # A caller whose step failed says so, not 'not run'.
    failed = {k: f"failed (step {sections[k]['step']})" for k in ("cpic", "pypgx", "cyrius")
              if sections[k]["state"] == "failed"}
    pg_ok = pg is not None and pgx_outside_calls.is_call(pg)
    cy_ok = cy is not None and pgx_outside_calls.is_call(cy.get("genotype")) and cy.get("filter") == "PASS"
    cy_txt = None
    if cy is not None:
        cy_txt = cy.get("genotype") or "none"
        if pgx_outside_calls.is_call(cy_txt) and not cy_ok:
            cy_txt += f" (not usable: Filter {cy.get('filter') or 'missing'})"
    calls = {
        "PharmCAT": (pc["diplotype"] if pc["status"] not in ("not called", "ambiguous") else pc["status"]) if pc
        else failed.get("cpic"),
        "pypgx": failed.get("pypgx", pg),
        "Cyrius": failed.get("cyrius", cy_txt),
    }
    agree = pgx_outside_calls.same_diplotype(pg, cy["genotype"]) if pg_ok and cy_ok else None
    try:
        consensus = read_consensus(sample_dir, sample, sections)
    except (OSError, csv.Error) as e:
        consensus = {"result": "unreadable", "passed_to_pharmcat": False, "reason": str(e), "source": None}
    # pypgx's call alone, labelled, when step 36's rule allows it
    # (pgx_outside_calls.probable_cyp2d6); it never reaches PharmCAT.
    probable = None
    if pg_ok:
        pd = sections["pypgx"]["data"]
        ev = pd.get("cyp2d6_evidence") or {}
        probable = pgx_outside_calls.probable_cyp2d6(
            {"Diplotype": pg, "Phenotype": ev.get("phenotype"), "Defining_filter": ev.get("defining_filter"),
             "Quality": ev.get("quality")},
            (cy.get("genotype"), cy.get("filter")) if cy is not None else None,
            (pd.get("depth_check") or {}).get("status"))
    return {"calls": calls, "agree": agree, "consensus": consensus, "probable": probable}


def collect(sample, sample_dir, declared_sex=None):
    status = read_run_status(sample_dir)
    manifest = read_manifest(sample_dir)
    sections = {}
    for key, title, step, reader in SECTIONS:
        sec = {"title": title, "step": step, "state": "missing", "source": None,
               "file_date": None, "note": None, "data": {}}
        try:
            path, data = reader(sample_dir, sample)
        except Unreadable as e:
            sec["state"], sec["note"] = "unreadable", str(e)
            sections[key] = sec
            continue
        except (OSError, ValueError, KeyError, IndexError, EOFError, csv.Error, zlib.error) as e:
            sec["state"], sec["note"] = "unreadable", f"{type(e).__name__}: {e}"
            sections[key] = sec
            continue
        if not path and status and (status["steps"].get(step) or "").startswith("failed"):
            # run-all.sh records 'failed' for a step it ran that exited non-zero
            sec["state"] = "failed"
            sec["note"] = (f"step {step} failed in the run of {status['started_utc'] or 'unknown date'}; "
                           "its log is in logs/")
        if path:
            mtime = os.path.getmtime(path)
            sec.update(state="ok", source=os.path.relpath(path, sample_dir), file_date=iso(mtime), data=data)
            if status and status["started_epoch"] is not None and step in RUN_ALL_STEPS:
                result = status["steps"].get(step)
                if result != "ok" and mtime < status["started_epoch"]:
                    sec["state"] = "stale"
                    sec["note"] = (f"from an earlier run ({iso(mtime)}); step {step} was "
                                   f"{result or 'not run'} in the run of {status['started_utc']}")
        sections[key] = sec

    cyp2d6 = cyp2d6_block(sample, sample_dir, sections)
    # The repeat readings lean on the traits (R117H) and CPIC (UGT1A1)
    # tables, whose sections come after the expansions one: read them again
    # with those sections' states.
    ex = sections["expansions"]
    if ex["state"] in ("ok", "stale"):
        ex["data"]["readings"] = str_readings({l["locus"]: l for l in ex["data"]["loci"]}, sample_dir, sample,
                                              {k: sections[k]["state"] for k in ("traits", "cpic")})
    if sections["locityper"]["state"] in ("ok", "stale"):
        primary = locityper_primary(sections)
        for x in sections["locityper"]["data"]["loci"]:
            x["primary"] = primary.get(x["locus"], "no other step calls this locus")
    declared = (declared_sex or (status or {}).get("declared_sex") or manifest_value(manifest, "run", "declared_sex")
                or sections["sample_qc"]["data"].get("declared_sex") or None)
    cnvpytor_sex_filter(sections["cnvpytor"]["data"], declared)

    clinvar_date = manifest_value(manifest, "data", "clinvar_file_date")
    summary = {
        "schema_version": SCHEMA_VERSION,
        "sample": sample,
        "generated_utc": iso(datetime.now(tz=timezone.utc).timestamp()),
        "run": {
            "status_file": status is not None,
            "started_utc": status["started_utc"] if status else None,
            "declared_sex": declared,
        },
        "databases": {
            "clinvar_release": clinvar_date,
            "cpsr_clinvar_release": manifest_value(manifest, "data", "cpsr_clinvar_release"),
            "hla_database": manifest_value(manifest, "data", "hla_database"),
        },
        "manifest": manifest or [],
        "sections": sections,
        "cyp2d6": cyp2d6,
        # The read count is written by alignment, which a run from a BAM or
        # a VCF never does, so it is no step left undone.
        "not_run": [sections[k]["title"] + (f" (step {sections[k]['step']} failed)" if sections[k]["state"] == "failed" else "")
                    for k, *_ in SECTIONS if sections[k]["state"] in ("missing", "failed") and k not in NOT_A_STEP],
        "not_assessed": NOT_ASSESSED + ([] if sections["prs"]["data"].get("adjusted") else [PRS_NOT_ADJUSTED])
        # Step 35b on a region subset (--paraphase_regions) may leave hba out
        + ([] if {"HBA1", "HBA2"} <= {g["gene"] for g in sections["paralogs"]["data"].get("genes", [])
                                      if g["region"] == "hba" and g["status"] == "called"}
           or (sections["parascopy"]["data"].get("hba") or {}).get("status") == "called"
           else [HBA_NOT_ASSESSED])
        + ([MEI_NOT_ASSESSED] if sections["mei"]["state"] == "missing" else [])
        + ([] if sections["vntyper"]["state"] == "ok"
           and sections["vntyper"]["data"].get("status") in ("positive", "possible", "negative")
           else [VNTYPER_NOT_ASSESSED])
        + ([] if sections["mosaic"]["data"].get("status") == "assessed" else [MCA_NOT_ASSESSED]),
    }
    recheck = run_delta.read_recheck(os.path.join(sample_dir, "recheck", "state.json"))
    if recheck:
        summary["recheck"] = recheck
    since = new_since(sample_dir, summary)
    if since:
        summary["new_since"] = since
    return summary


def new_since(sample_dir, summary):
    """{previous_utc, sections}: what changed since the newest run kept in history/ (generate-report.sh
    keeps one per run), or None without one. A folder that cannot be read counts as none."""
    prev_dir = run_delta.newest_history(sample_dir)
    if not prev_dir:
        return None
    try:
        with open(os.path.join(prev_dir, "summary.json")) as f:
            prev = json.load(f)
        changes = run_delta.compare(prev, run_delta.read_manifest(os.path.join(prev_dir, "run_manifest.tsv")),
                                    summary, run_delta.read_manifest(os.path.join(sample_dir, "run_manifest.tsv")))
    except (OSError, ValueError, AttributeError, TypeError, KeyError):
        return None
    return {"previous_utc": prev.get("generated_utc"), "sections": run_delta.changed_sections(changes)}


def sample_qc_main(argv):
    ap = argparse.ArgumentParser(prog="collect_summary.py sample-qc",
                                 description="Write the step 33 table (key, value) for one sample.")
    ap.add_argument("--sample", required=True)
    ap.add_argument("--somalier-samples", required=True, help="somalier relate's samples.tsv")
    ap.add_argument("--somalier-pairs", help="somalier relate's pairs.tsv (other samples of the run)")
    ap.add_argument("--somalier-id", help="the sample's name in somalier's files (the BAM's @RG SM in the bash step, the samplesheet id in Nextflow; default --sample)")
    ap.add_argument("--selfsm", help="VerifyBamID2's .selfSM")
    ap.add_argument("--declared-sex", choices=["male", "female", ""], default="")
    ap.add_argument("--freemix-warn", type=float, default=FREEMIX_WARN)
    ap.add_argument("--marker-check", choices=["passed", "skipped", ""], default="",
                    help="skipped: VerifyBamID2 ran with --DisableSanityCheck (fewer than 1,000 markers had reads)")
    ap.add_argument("--out", required=True)
    a = ap.parse_args(argv)
    try:
        rows = sample_qc_table(a.sample, a.somalier_samples, a.selfsm, a.somalier_pairs, a.declared_sex,
                               a.freemix_warn, a.somalier_id, a.marker_check)
    except (Unreadable, OSError) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        return 1
    tmp = a.out + ".tmp"
    with open(tmp, "w") as f:
        f.write("key\tvalue\n")
        for k, v in rows:
            f.write(f"{k}\t{v}\n")
    os.replace(tmp, a.out)
    for k, v in rows:
        print(f"  {k}: {v}")
    return 0


def main(argv=None):
    argv = sys.argv[1:] if argv is None else argv
    if argv[:1] == ["sample-qc"]:
        return sample_qc_main(argv[1:])
    if argv[:1] == ["prs-format"]:
        return prs_format_main(argv[1:])
    if argv[:1] == ["prs-table"]:
        return prs_table_main(argv[1:])
    if argv[:1] == ["mito-annotate"]:
        return mito_annotate_main(argv[1:])
    if argv[:1] == ["mito-deletions"]:
        return mito_deletions_main(argv[1:])
    if argv[:1] == ["mito-copies"]:
        return mito_copies_main(argv[1:])
    ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
    ap.add_argument("--sample", required=True)
    ap.add_argument("--sample-dir", required=True)
    ap.add_argument("--declared-sex", help="the samplesheet's sex, when no run status or manifest records it")
    ap.add_argument("--out", help="summary JSON (default: SAMPLE_DIR/summary.json)")
    a = ap.parse_args(argv)
    summary = collect(a.sample, a.sample_dir, a.declared_sex)
    out = a.out or os.path.join(a.sample_dir, "summary.json")
    tmp = out + ".tmp"
    with open(tmp, "w") as f:
        json.dump(summary, f, indent=1)
        f.write("\n")
    os.replace(tmp, out)
    print(f"Summary written: {out}")
    return 0


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