💡 導讀：地圖上的第 100 號門牌，為什麼有時寫成 99？

Day 17，我們學會從歷史 AlphaGenome Atlas AVI 產物讀取單點變異線索；Day 18，從歷史 Model REF／ALT 輸出認識剪接差異；Day 19，以五筆人工查詢練習 asyncio 與錯誤重試。Day 20 的問題更基本：當我們看到變異 chr1:111:A>G，要如何判斷它落在哪個研究用的區間標籤？

就像兩種地圖可能把同一條街的門牌從 0 或 1 開始編號，VCF 與 BED 也有兩種坐標制度。如果沒先換算，分析可能錯一個鹼基；更糟的是，把「位置重疊」直接當成「變異破壞轉錄因子結合」。今天先把這兩件事分清楚，再學習怎麼從 AI Studio 面板逆向拆解成 Python。

一、核心概念：VCF、BED 與半開區間

換算規則是：1-based 封閉區間 [start,end] → 0-based 半開區間 [start−1,end)。以人工區間 A chr1:[99,104)、B chr1:[109,114)、C chr1:[120,126)、D chr2:[99,103) 為例：

判斷重疊時要先確認染色體一致，再使用 max(0, min(end1,end2) - max(start1,start2))。重疊只是資料定位，不能證明 motif 序列確實存在、蛋白質已結合、變異造成致病。 這四個 A/B/C/D 標籤是人工設計，並非真的執行了 motif 序列掃描。

二、Google AI Studio：先做出 Motif 區間查詢面板

程式小白可以先把下面的提示詞貼到 Google AI Studio Build Mode，做出純展示用的研究 MVP，再對照今天的 Python 函式名稱。MVP 使用人工資料，不需要或不應要求使用者輸入 API Key。

請用繁體中文建立「基因座區間探索」Research Use Only 示範面板。
所有資料都是 SYNTHETIC EDUCATIONAL DEMO，不連任何後端服務。

輸入區：
- 1-based 的區間輸入框，預設 chr1:100-120
- 教學用單點 SNV 輸入框，預設 chr1:111:A>G
- 四個 BED 樣式人工標籤 A、B、C、D（全部標為 SYNTHETIC）

輸出區：
- 清楚顯示 1-based closed → 0-based half-open 的座標換算
- A、B 命中；C 因右端點不包含而排除；D 因染色體不同而排除
- 高亮人工 SNV 111 只落在 B
- 清楚顯示 overlap_bp 和查詢結果 JSON

禁止輸出臨床疾病風險、真實 AVI 分數、官方 Motif Atlas 結果。
禁止外部 API、病人資料、金鑰與登入操作；標示「區間重疊 ≠ TF binding」。

畫出來後，再用「逆向拆解學習法」檢查：

AI Studio 產生的介面通常仍需你自己理解資料流與後端呼叫；不是按 Get Code 就能自動得到這篇 Python 程式。這篇實作的是與介面對應的離線資料處理核心。

💻 三、今天動手做：獨立可跑的 Python Mini Demo

請將下列完整程式存成 day20_offline_motif_lab.py，並將本教學包的 day20_synthetic_motifs.tsv 放在同一資料夾。只需要 Python 3.12 標準函式庫。程式只允許 --mode offline，即使電腦上已設定 API Key，也不會呼叫任何網路服務。

"""Day 20: OFFLINE-only synthetic motif/interval coordinate learning lab.

All motif labels are invented BED-like intervals, not motif sequence hits.
No AlphaGenome SDK, network access, patient data, or real API calls.
Python 3.12+ standard library only.
"""
from __future__ import annotations

import argparse
import csv
from dataclasses import asdict, dataclass
from html import escape
import json
from pathlib import Path
import re

REGION_PATTERN = re.compile(r"^(chr(?:[1-9]|1[0-9]|2[0-2]|X|Y|M|Toy)):(\d+)-(\d+)$")
VARIANT_PATTERN = re.compile(
    r"^(chr(?:[1-9]|1[0-9]|2[0-2]|X|Y|M|Toy)):(\d+):([ACGT])>([ACGT])$"
)
MOTIF_ID_PATTERN = re.compile(r"^[A-Za-z0-9][A-Za-z0-9._-]{0,63}$")
MAX_WINDOW_BP = 10_000


@dataclass(frozen=True)
class MotifInterval:
    motif_id: str
    chrom: str
    start_0based: int
    end_0based: int
    evidence_class: str = "SYNTHETIC_FIXTURE"


DEFAULT_MOTIFS = (
    MotifInterval("A", "chr1", 99, 104),
    MotifInterval("B", "chr1", 109, 114),
    MotifInterval("C", "chr1", 120, 126),
    MotifInterval("D", "chr2", 99, 103),
)


def valid_motifs(rows):
    """Reject invalid or non-synthetic intervals instead of silently accepting them."""
    rows = tuple(rows)
    if not rows or len(rows) > 1000:
        raise ValueError("Require 1..1000 synthetic motif intervals")
    ids = set()
    for row in rows:
        if not isinstance(row, MotifInterval):
            raise ValueError("Motif row type must be MotifInterval")
        if not MOTIF_ID_PATTERN.fullmatch(row.motif_id):
            raise ValueError("Motif ID format is invalid")
        if row.motif_id in ids:
            raise ValueError("Duplicate motif ID")
        ids.add(row.motif_id)
        if not REGION_PATTERN.fullmatch(f"{row.chrom}:1-1"):
            raise ValueError("Motif chromosome format is invalid")
        if (type(row.start_0based) is not int or type(row.end_0based) is not int
                or row.start_0based < 0 or row.end_0based <= row.start_0based):
            raise ValueError("Motifs must use nonempty 0-based half-open intervals")
        if row.evidence_class != "SYNTHETIC_FIXTURE":
            raise ValueError("Only expressly synthetic fixtures are accepted")
    return rows


def parse_region_1based(value):
    """1-based closed region string -> chromosome and 0-based [start,end)."""
    match = REGION_PATTERN.fullmatch(value)
    if not match:
        raise ValueError("Region must have form chr1:100-120 (1-based, closed)")
    chrom, raw_start, raw_end = match.groups()
    start, end = int(raw_start), int(raw_end)
    if start < 1 or end < start or (end - start + 1) > MAX_WINDOW_BP:
        raise ValueError("Region must be positive, ordered, and <= 10000 bases")
    return chrom, start - 1, end


def parse_variant_1based(value):
    match = VARIANT_PATTERN.fullmatch(value)
    if not match:
        raise ValueError("Teaching SNV must look like chr1:111:A>G")
    chrom, pos, ref, alt = match.groups()
    position = int(pos)
    if position < 1 or ref == alt:
        raise ValueError("Teaching SNV requires position >=1 and REF != ALT")
    return chrom, position - 1, ref, alt


def overlap_bp(a_start, a_end, b_start, b_end):
    return max(0, min(a_end, b_end) - max(a_start, b_start))


def load_synthetic_tsv(path):
    """Load user-supplied explicit synthetic fixture, without any genomic lookup."""
    with Path(path).open("r", encoding="utf-8-sig", newline="") as handle:
        reader = csv.DictReader(handle, delimiter="\t")
        needed = ["motif_id", "chrom", "start_0based", "end_0based", "evidence_class"]
        if reader.fieldnames != needed:
            raise ValueError("Fixture TSV columns must exactly match the documented header")
        out = []
        for row in reader:
            try:
                out.append(MotifInterval(
                    row["motif_id"], row["chrom"], int(row["start_0based"]),
                    int(row["end_0based"]), row["evidence_class"]
                ))
            except (TypeError, ValueError) as exc:
                raise ValueError("Fixture coordinates must be integers") from exc
    return valid_motifs(out)


def analyze(region="chr1:100-120", variant="chr1:111:A>G", motifs=DEFAULT_MOTIFS):
    motifs = valid_motifs(motifs)
    chrom, query_start, query_end = parse_region_1based(region)
    variant_chrom, variant_zero, ref, alt = parse_variant_1based(variant) if variant is not None else (None, None, None, None)
    variant_in_query = bool(variant_chrom == chrom and query_start <= variant_zero < query_end)
    overlaps = []
    for item in motifs:
        if item.chrom != chrom:
            continue
        n = overlap_bp(query_start, query_end, item.start_0based, item.end_0based)
        if n:
            overlaps.append({**asdict(item), "overlap_bp": n,
                             "contains_variant": bool(variant_in_query and item.start_0based <= variant_zero < item.end_0based)})
    return {
        "evidence_class": "SYNTHETIC_EDUCATIONAL_DEMO",
        "mode": "OFFLINE_ONLY",
        "query_1based_closed": region,
        "query_0based_half_open": {"chrom": chrom, "start": query_start, "end": query_end},
        "variant_1based": variant,
        "variant_in_query": variant_in_query,
        "variant_overlap_ids": [r["motif_id"] for r in overlaps if r["contains_variant"]],
        "overlap_count": len(overlaps),
        "overlap_ids": [r["motif_id"] for r in overlaps],
        "region_overlaps": overlaps,
        "motif_interval_count": len(motifs),
        "motif_sequence_scanned": False,
        "reference_allele_verified": False,
        "official_motif_atlas": False,
        "atlas_query_performed": False,
        "clinical_interpretation": False,
        "real_api_calls_this_run": 0,
    }


def render_html(result):
    """Produce a static offline teaching panel with escaped user-facing strings."""
    rows = []
    for item in result["region_overlaps"]:
        cells = (item["motif_id"], item["chrom"],
                 f'[{item["start_0based"]}, {item["end_0based"]})',
                 str(item["overlap_bp"]), "YES" if item["contains_variant"] else "NO")
        rows.append("<tr>" + "".join(f"<td>{escape(x)}</td>" for x in cells) + "</tr>")
    query = escape(result["query_1based_closed"])
    bed = result["query_0based_half_open"]
    bed_display = escape(f'{bed["chrom"]}:[{bed["start"]},{bed["end"]})')
    variant = escape(result["variant_1based"] or "not provided")
    table = "\n".join(rows) if rows else '<tr><td colspan="5">No overlaps in synthetic fixture</td></tr>'
    return f"""<!doctype html>
<html lang="zh-Hant"><head><meta charset="utf-8"><meta name="viewport" content="width=device-width,initial-scale=1">
<title>Day 20 | Offline motif intervals</title><style>
body{{font-family:system-ui,-apple-system,sans-serif;background:#E5DEDA;color:#27313b;margin:0;padding:2rem}}
main{{background:#fffdfb;border-radius:18px;max-width:900px;margin:auto;padding:2rem;box-shadow:0 8px 32px #0001}}
small{{color:#765f6b}} h1{{font-size:1.55rem}} .chip{{display:inline-block;background:#eaf4f7;padding:.35rem .7rem;border-radius:20px;margin:.2rem}}
table{{width:100%;border-collapse:collapse;margin-top:1.2rem}}th,td{{padding:.7rem;border-bottom:1px solid #ddd;text-align:left}}
th{{background:#f7f2f0}}.warn{{background:#fff2e8;padding:1rem;border-radius:10px;border-left:4px solid #bd7659}}
</style></head><body><main>
<small>SYNTHETIC EDUCATIONAL DEMO · OFFLINE_ONLY · no AlphaGenome API</small>
<h1>Day 20 · 基因座區間探索</h1>
<p><span class="chip">1-based: {query}</span> <span class="chip">0-based BED: {bed_display}</span></p>
<p>教學變異：<b>{variant}</b>；位於人工標籤中的候選：<b>{escape(", ".join(result["variant_overlap_ids"]) or "none")}</b></p>
<p>重疊人工區間共 <b>{result["overlap_count"]}</b> 筆。</p>
<table><thead><tr><th>Motif label</th><th>Chromosome</th><th>BED interval</th><th>Overlap (bp)</th><th>Contains SNV</th></tr></thead><tbody>{table}</tbody></table>
<p class="warn">這些是人工命名的 BED 標籤，不是 DNA motif 序列掃描結果，也不是蛋白質結合、AlphaGenome 官方 Atlas 資料或臨床判讀。</p>
</main></body></html>
"""


def run(output_dir="day20_output", region="chr1:100-120", variant="chr1:111:A>G", motifs=DEFAULT_MOTIFS):
    result = analyze(region=region, variant=variant, motifs=motifs)
    folder = Path(output_dir)
    folder.mkdir(parents=True, exist_ok=True)
    json_file = folder / "day20_interval_report.json"
    html_file = folder / "day20_motif_view.html"
    json_file.write_text(json.dumps(result, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
    html_file.write_text(render_html(result), encoding="utf-8")
    return {"status": "PASS", "mode": "OFFLINE_ONLY", "query": region,
            "overlap_count": result["overlap_count"], "overlap_ids": result["overlap_ids"],
            "variant_overlap_ids": result["variant_overlap_ids"],
            "json_file": str(json_file), "html_file": str(html_file),
            "real_api_calls_this_run": 0}


def main(argv=None):
    parser = argparse.ArgumentParser(description="Day 20 offline interval/motif label lesson")
    parser.add_argument("--mode", choices=["offline"], default="offline")
    parser.add_argument("--region", default="chr1:100-120", help="1-based inclusive query")
    parser.add_argument("--variant", default="chr1:111:A>G", help="Synthetic SNV; not REF-verified")
    parser.add_argument("--fixture", help="Optional SYNTHETIC_FIXTURE TSV")
    parser.add_argument("--output-dir", default="day20_output")
    args = parser.parse_args(argv)
    motifs = load_synthetic_tsv(args.fixture) if args.fixture else DEFAULT_MOTIFS
    print(json.dumps(run(args.output_dir, args.region, args.variant, motifs),
                     ensure_ascii=False, sort_keys=True))


if __name__ == "__main__":
    main()

▶ 執行方式（WSL Ubuntu）

如果程式下載到 Windows Downloads，在 WSL 輸入：

cd /mnt/c/Users/user/Downloads

python3.12 day20_offline_motif_lab.py \
  --mode offline \
  --region 'chr1:100-120' \
  --variant 'chr1:111:A>G' \
  --fixture day20_synthetic_motifs.tsv \
  --output-dir day20_output

程式會建立 day20_output/day20_interval_report.json 與 day20_output/day20_motif_view.html。雙擊 HTML 即可從瀏覽器查看人工 A/B 命中表，完全不需要雲端後端或網路連線。

▶ 預期與本機實測輸出

我在 Python 3.13.5 環境執行，Exit Code 0；標準輸出為：

{"html_file": "day20_output/day20_motif_view.html", "json_file": "day20_output/day20_interval_report.json", "mode": "OFFLINE_ONLY", "overlap_count": 2, "overlap_ids": ["A", "B"], "query": "chr1:100-120", "real_api_calls_this_run": 0, "status": "PASS", "variant_overlap_ids": ["B"]}

JSON 詳細結果中的重要部分：

{"query_0based_half_open": {"chrom":"chr1","start":99,"end":120}, "overlap_ids":["A","B"], "variant_overlap_ids":["B"], "motif_sequence_scanned":false, "official_motif_atlas":false, "atlas_query_performed":false, "real_api_calls_this_run":0}

它只證明人工 TSV 區間重疊：A、B 命中，變異位於 B。沒有真的對 chr1:111:A>G 核對 GRCh38 的 REF，也沒有呼叫模型預測。

🔎 右邊界小實驗

現在改查 chr1:121-121，將變異一起改成 chr1:121:A>G：

python3.12 day20_offline_motif_lab.py \
  --mode offline \
  --region 'chr1:121-121' \
  --variant 'chr1:121:A>G' \
  --fixture day20_synthetic_motifs.tsv \
  --output-dir day20_right_edge_output

我在 Python 3.13.5 實測 stdout：

{"html_file": "day20_right_edge_output/day20_motif_view.html", "json_file": "day20_right_edge_output/day20_interval_report.json", "mode": "OFFLINE_ONLY", "overlap_count": 1, "overlap_ids": ["C"], "query": "chr1:121-121", "real_api_calls_this_run": 0, "status": "PASS", "variant_overlap_ids": ["C"]}

由於位置 121 對應 BED [120,121)，它恰好命中 C，而不再命中 A／B。這個練習是 Day 20 最重要的 Checkpoint。

✅ 本日 Checklist

• [ ] 我能說明為什麼 1-based chr1:100-120 等於 BED [99,120)，長度都是 21 bp。
• [ ] 我能解釋 A／B 重疊、C 因端點不包含而排除、D 因染色體不同而排除。
• [ ] 我知道 SNV chr1:111:A>G 在人工表裡落入 B，只是座標標註而不是生物機制。
• [ ] 我能開啟本機 HTML 並查看 JSON 裡 overlap_ids、variant_overlap_ids 與 real_api_calls_this_run=0。
• [ ] 我能把查詢改成 chr1:121-121，得到 ['C']。
• [ ] 我能分清楚「本地 Motif 區間清單」、「官方 Atlas 預計算 SNV 分數」與「真正 TF binding／motif 作用實驗」。

四、向真正 Atlas 區間查詢銜接：只先看懂 API 介面

Google DeepMind 的 Atlas API 有 query_interval()，其輸入使用 genome.Interval，座標採 0-based、右端不包含。當得到正式存取權、確認 SDK 版本，且作者另外批准 LIVE 查詢後，才可查指定區間內的 AVI 等預計算分數。該 API 不會因我們今天的 A/B/C/D TSV 標籤而自動證明 motif 結合。

以下只是文件介面片段，沒有放入可執行的教學程式，也未執行：

# OFFICIAL INTERFACE REFERENCE ONLY — NOT EXECUTED
# After separate authorization, a verified AtlasClient could call:
# results = client.query_interval(
#     genome.Interval(chromosome="chr1", start=99, end=120),
#     requested_scorers=["AVI_SCORE"],
# )
# This queries precomputed SNV impact scores, NOT motif hits.

Day 17 的 Reboot 有 SDK 0.9.0、query_variant() 歷史成功證據；本教學包未測 query_interval() LIVE。Bridge 0.10.0 的 SDK 安裝版／LIVE 相容性也尚未完成驗證。不能將 Day 17 歷史單點成功記錄說成 Day 20 區間 API 已成功。所有正式輸出還需核對 genome build、參考序列、計分方法及授權條款。

官方技術參考：

• AlphaGenome Atlas API 官方索引
• Google DeepMind 官方 Atlas 變異評分開發文件
• AlphaGenome 基因組座標資料型別

五、科學證據與限制

SYNTHETIC EDUCATIONAL DEMO：本篇人工 A/B/C/D 並非任何正式 PWM、Motif Atlas、ATAC、ChIP-seq 或 TF 結合實驗。它沒有序列掃描，不能判斷位點是否含特定 motif，也不能推論某個變異會造成蛋白質結合改變、致病或病患風險。人工 chr1 座標也不是經 GRCh38 FASTA 驗證的真實 SNV。

OFFLINE SOFTWARE QA：新的腳本在 Python 3.13.5 執行，22/22 單元測試 PASS；本次工作環境沒有 Python 3.12，因此不能冒稱已在此環境完成 Python 3.12 實測。請讀者在自己的 WSL 3.12 環境重跑。

CURRENT LIVE API = NOT_RUN：沒有 AlphaGenome／Gemini／PubMed API 請求，沒有金鑰，沒有查詢臨床個案；真實 API 需要獨立的存取授權與版本確認。GitHub 保持 PRIVATE，PUBLICATION_ON_HOLD。

六、明日預告

明天 Day 21，我們將從 VCF 到研究審閱封包：學習讓變異經過 QC、註釋與證據整理，再交由研究者審閱，而不是僅因某個位置命中標籤就判定具有疾病意義。