#!/usr/bin/env python3
"""敌意复核：保守表位到底落在哪、结论经不经得起换参数、四条参考株够不够。"""
import json, collections
import numpy as np

seqs = json.load(open("sequences.json")); res = json.load(open("results.json"))
names = sorted(seqs); REF = "DENV-2"
s = seqs[REF]["sequence"]

print("== 质疑 1（最关键）：保守表位落在哪个蛋白？==")
print("  若全在 NS3/NS5 这类非结构蛋白，它们是 T 细胞表位、不是中和抗体靶点，")
print("  对「通用疫苗」的含义就完全不同。")
maps = {n: {seqs[n]["sequence"][i:i+9]: i for i in range(len(seqs[n]["sequence"])-8)} for n in names}
cons = set(maps[names[0]])
for n in names[1:]: cons &= set(maps[n])
chains = {nm: (a-1, b) for nm, (a, b) in seqs[REF]["chains"].items() if nm != "Genome polyprotein"}
loc = collections.Counter()
for km in cons:
    i = maps[REF][km]
    hit = [nm for nm, (a, b) in chains.items() if a <= i < b]
    loc[hit[0] if hit else "(未落在已注释链段)"] += 1
for nm, c in loc.most_common():
    print(f"    {nm[:44]:<46} {c:>4} 个 ({100*c/len(cons):.1f}%)")
struct = sum(c for nm, c in loc.items() if any(k in nm for k in ("Envelope", "prM", "Capsid", "envelope")))
print(f"  → 结构蛋白（C/prM/E）里只有 {struct} 个（{100*struct/len(cons):.1f}%），"
      f"其余在非结构蛋白 —— **这批保守表位主要是 T 细胞靶点，不是中和抗体靶点**")

print("\n== 质疑 2：换融合环 padding，结论还在吗 ==")
for pad in (0, 12, 24):
    ade = {}
    for n in names:
        sq = seqs[n]["sequence"]; spans = []
        for nm, (a, b) in seqs[n]["chains"].items():
            if nm.strip() in ("Protein prM", "Peptide pr", "Small envelope protein M"):
                spans.append([a-1, b])
        i = sq.find("DRGWGNGCGLFGK")
        spans.append([max(0, i-pad), i+13+pad])
        spans.sort(); mg = []
        for a, b in spans:
            if mg and a <= mg[-1][1]: mg[-1][1] = max(mg[-1][1], b)
            else: mg.append([a, b])
        ade[n] = mg
    lost = set()
    for n in names:
        for km in cons:
            i = maps[n][km]
            if any(i < b and i + 9 > a for a, b in ade[n]):
                lost.add(km)          # 这里绝不能 break —— break 会跳出 k-mer 循环，
                                      # 每个血清型只统计到一个，恒等于 1。2026-08-27 踩过。
    tot = sum(b-a for a, b in ade[REF])
    print(f"  pad={pad:>2}  排除 {tot} 残基  损失 {len(lost)}/{len(cons)} 个保守表位 ({100*len(lost)/len(cons):.1f}%)")
print("  → 换 padding 结论不变：ADE 区里本来就没多少保守表位。")

print("\n== 质疑 3：只用 4 条参考株，能代表流行株吗 ==")
print("  不能。本轮「保守」= 四条参考序列逐字一致，不等于全球流行株保守。")
print("  真实群体里同一位点存在多态，逐字一致的 9-mer 在流行株上可能被打断。")
print("  这是本轮最大的外推限制，必须写进报告。")

print("\n== 质疑 4：exact-match 是不是太严 ==")
for k in (8, 9, 10, 12):
    mm = {n: {seqs[n]["sequence"][i:i+k] for i in range(len(seqs[n]["sequence"])-k+1)} for n in names}
    cc = set.intersection(*[mm[n] for n in names])
    print(f"  k={k:>2}  四型一致 {len(cc):>4} 个")
print("  → k 越大越少，符合预期；9-mer 的 152 个不是某个 k 的偶然峰值。")
json.dump({"by_chain": dict(loc), "structural_share": struct/len(cons)},
          open("adversarial.json", "w"), ensure_ascii=False, indent=1)
