# Transition share and hotspot share of PhyloTree Build 17 substitutions, binned by depth below the root.
# Run from research/:  python3 code/depth_spectrum.py code data   (copied from a scratch analysis on 2026-10-27;  port to Deno pending).
import sys, re, csv, os
CODE = sys.argv[1]; DATA = sys.argv[2]
sys.path.insert(0, CODE)
import phylotree17 as pt
root, idx, allnodes = pt.load()
rec = {}
with open(os.path.join(DATA, 'phylotree-b17-position-recurrence.tsv')) as f:
    for r in csv.DictReader(f, delimiter='\t'):
        rec[int(r['rcrs_position'])] = int(r['n_independent_branches'])
def toks(muts):
    out = []
    for raw in muts.split():
        if '(' in raw or ')' in raw: continue
        t = raw.rstrip('!')
        if t.endswith('d') or '.' in t or re.match(r'^\d+-\d+', t): continue
        m = re.match(r'^([ACGT])?(\d+)([ACGTacgt])$', t)
        if not m: continue
        pos = int(m.group(2)); der = m.group(3)
        out.append((pos, der.islower()))
    return out
bins = {}
def start_depth(n):
    return pt.depth_subs(n['parent']) if n['parent'] is not None and n['parent']['depth'] >= 0 else 0
for n in pt.walk(root):
    if n['depth'] < 0: continue
    sd = start_depth(n)
    b = '0-9' if sd < 10 else '10-19' if sd < 20 else '20-29' if sd < 30 else '30-39' if sd < 40 else '40+'
    leaf = 'tip' if not n['kids'] else 'internal'
    for key in (b, 'ALL', leaf):
        d = bins.setdefault(key, [0, 0, 0, 0, 0])
        for pos, tv in toks(n['muts']):
            d[0] += 1; d[1] += tv
            r = rec.get(pos, 0)
            d[2] += (r >= 21); d[3] += (r <= 2)
            d[4] += (pos <= 576 or pos >= 16024)
print('bin\tsubs\ttransversion%\thotspot(>=21)%\trare(<=2)%\tcontrol%')
for k in ['0-9','10-19','20-29','30-39','40+','internal','tip','ALL']:
    d = bins.get(k)
    if not d: continue
    print(f"{k}\t{d[0]}\t{100*d[1]/d[0]:.1f}\t{100*d[2]/d[0]:.1f}\t{100*d[3]/d[0]:.1f}\t{100*d[4]/d[0]:.1f}")
# control-region-only and coding-only transversion share, deep (<30 subs from root) vs shallow
agg = {}
for n in pt.walk(root):
    if n['depth'] < 0: continue
    sd = start_depth(n)
    zone = 'deep<30' if sd < 30 else 'shallow>=30'
    for pos, tv in toks(n['muts']):
        reg = 'control' if (pos <= 576 or pos >= 16024) else 'coding'
        d = agg.setdefault((zone, reg), [0, 0]); d[0] += 1; d[1] += tv
for k in sorted(agg):
    d = agg[k]; print(k, d[0], f"{100*d[1]/d[0]:.1f}% transversions")