系统发育进化树构建

> 把「拿一批序列 → 出可投稿的进化树」从反复排障的一天，变成改一个 RECORDS
> 列表加一条命令。

从一组序列（或基因名/accession）出发，产出带 bootstrap 支持度的进化树图。
另含两条基因家族分析常用扩展：基因结构图、远缘外群定根。

什么时候用

• 「对 X 基因/蛋白做系统发育分析」「建进化树」「看看这几个基因的亲缘关系」
• 需要判断一组基因是直系同源（ortholog）还是旁系同源（paralog）
• 需要给 Gene family 分析出 phylogenetic tree 图（论文 Figure）
• 手上有一批 accession / FASTA，要出树
• 「把序列换成基因结构图」「标注外显子/内含子」「加上蛋白结构域」→ 工作流 A
• 「加外群定根」「判断哪个先分化」「无根树要定根」→ 工作流 B

什么时候不用

• 只想要单条序列的注释（结构域/理化性质）→ 用 protein-seq-annotation
• 只要引物设计 → 用 gene-primer-design
• 要做全基因组比较基因组学（共线性/Ka-Ks）→ 本 skill 只做树，不含这些
• 只需要 qPCR 的 ΔΔCt 计算 → 用 qpcr-ddct-analysis

核心原则

1. 数字必须可追溯。 树上的支持度、分支长度、同一性都要能说明是怎么算的：
用的什么距离模型、什么过滤阈值、多少次 bootstrap。写进图注和报告。

2. 自写建树必须交叉验证。 这是本 skill 最重要的一条。自实现 NJ/UPGMA 最阴的
失败模式是「拓扑看着对、分支归属错」或「序列被静默混同」——
肉眼看图完全发现不了。必须做两道检查：

• 对齐后去 gap 的序列，必须与原序列逐字一致（断言）
• 自实现拓扑 vs Biopython 官方 DistanceTreeConstructor().nj()，
  参考实现的每个分支都必须被复现（NJ 在 Q 值并列时允许我方多解出 1-2 个支）

3. 用 accession 而不是基因名取序列。 基因名在不同文献里会漂移
（水稻 CAT 家族就同时存在 CATA/B/C 与 OsCAT1/2/3 两套编号，且对应关系不一致）。
按 accession 取，结果才是确定的；输出里同时打 accession + 基因座号，方便对文献。

4. 不要编拓扑。 树跑出来什么样就报什么样。如果外类群没形成干净的单系，
如实在报告里说明（常见原因：谱系分选不完全 ILS、外群选得太近），
不要为了「好看」调整或隐瞒。

快速开始

# 1) 装依赖到隔离环境（一次性）
PY="C:/Users/<user>/.workbuddy/binaries/python/versions/3.13.12/python.exe"
"$PY" -m venv "C:/Users/<user>/.workbuddy/binaries/python/envs/default"
VENV="C:/Users/<user>/.workbuddy/binaries/python/envs/default"
"$VENV/Scripts/pip.exe" install -q biopython numpy matplotlib

# 2) 把模板脚本复制到项目并改名
mkdir -p myphylo && cp references/phylo_pipeline.py myphylo/my_tree.py
cd myphylo

# 3) 改 my_tree.py 顶部的 RECORDS，然后跑
"$VENV/Scripts/python.exe" my_tree.py --offline --bootstrap 1000

期望输出（节选，实测）：

[联网] 拉取 16 条序列 ...
[2/5] 多序列比对  -> 493 列
[3/5] 位点过滤（非 gap >= 0.9）-> 492 列
[4/5] 距离矩阵（p-distance）
      OsCATA vs OsCATB = 72.15%   OsCATA vs OsCATC = 75.41%
      OsCATB vs OsCATC = 83.94%
[5/5] NJ + UPGMA 建树（bootstrap 1000x）
[OK] figures/my_tree_NJ.pdf  figures/my_tree_NJ.png
[OK] figures/my_tree_UPGMA.pdf  figures/my_tree_UPGMA.png

前置条件

环境准备（隔离，别污染用户环境）

PY="C:/Users/<user>/.workbuddy/binaries/python/versions/3.13.12/python.exe"
"$PY" -m venv "C:/Users/<user>/.workbuddy/binaries/python/envs/default"
VENV="C:/Users/<user>/.workbuddy/binaries/python/envs/default"
"$VENV/Scripts/pip.exe" install -q biopython numpy matplotlib

验证："$VENV/Scripts/python.exe" -c "import Bio,numpy,matplotlib;print(Bio.__version__)"

工作流

1. 明确输入与范围

先问清楚（或从上下文确认）：
- 哪些序列/基因？（要 accession 或 FASTA）
- 要不要参考同源基因？要不要外类群？
- 只跑目标基因（3-5 条）→ 树没定根方向，只能看相对距离
- 加同科/同属同源基因 → 能判断直系/旁系同源（推荐）
- 加远缘外群（如苔藓、蕨类、酵母）→ 才能定根
- 出图要不要投期刊？（要 → 出 PDF，且设 pdf.fonttype=42）

2. 取序列

UniProt REST（按 accession，最可靠）：

Q="accession:Q0E4K1+OR+accession:Q0D9C4+OR+accession:Q10S82"
curl -s "https://rest.uniprot.org/uniprotkb/stream?query=$Q&format=fasta" -o refs.fasta

按基因名/物种搜（用于找同源基因）：

curl -s "https://rest.uniprot.org/uniprotkb/search?\
query=protein_name:%22Catalase%22+AND+organism_id:3702+AND+reviewed:true\
&format=tsv&fields=accession,id,gene_names,length&size=20"

常用 organism_id：水稻(日本晴) 39947、拟南芥 3702、玉米 4577、大麦 4513、大豆 3847。

务必核对：长度是否合理（同一家族内不应差太多）、有没有注释成别的蛋白
（基因名会在不同物种间撞车，如 CAT1 在人类是阳离子氨基酸转运蛋白）。

3. 比对 → 建树 → 出图

直接把 references/phylo_pipeline.py 复制到项目里当脚本用，它已包含全部护栏。
只需改文件顶部的 RECORDS 列表：

RECORDS = [
    dict(acc="Q0E4K1", label="OsCATA", group="OsCAT",
         sp="Oryza sativa", locus="LOC_Os02g02400", note="Catalase A"),
    ...
]

跑（复制到项目后按你要的名字调用，下例叫 my_tree.py）：

"$VENV/Scripts/python.exe" my_tree.py --offline --bootstrap 1000

references/phylo_pipeline.py 已实现：取序列（UniProt，带本地缓存）、星型渐进式
MSA、保守位点过滤、p-distance 距离矩阵、NJ + UPGMA、bootstrap、PDF/PNG 出图，
以及全部自检断言。不要重写它 —— 直接改 RECORDS 和输出前缀。

4. 出图后一定要看图

用 Read 工具打开 PNG 亲眼确认：
- 叶标签有没有和分支线重叠 → 调 label_x 留白
- 支持度有没有超出 0-100%（超出 = 百分比换算 bug）
- 拓扑是否符合生物学预期
- 图例有没有遮挡标签

这一步不能省——脚本跑通 ≠ 图能用。

5. 报告结论

给出：
- 两两同一性表（目标基因之间）
- 主要分支 + bootstrap 支持度
- 生物学解读（直系/旁系同源？复制事件几次？）
- 不确定性说明（哪些分支支持度低、外群是否干净）

参数选择速查

Biopython 模型名坑：p-distance 在 Biopython 里叫 identity；
DistanceCalculator 支持的模型是 identity/blastp/blosum62/tanaka/kimura/pam250/jones/dayhoff/...
——没有 p-distance 这个名字，直接传会 ValueError。

已知坑（都已在模板脚本里加护栏）

见 references/pitfalls.md，包含每个坑的症状、根因与修法。三条最要命的：

1. MSA 合并时中心序列少算插入列 → 去 gap 后长度对不上。
2. 插入块顺序放反 → 长度对、内容错（更隐蔽）。
3. bootstrap 百分比重复乘 100 → 出现 10000% 的支持度。

外加两个作用域坑：--core-only 时筛选后的列表没真正生效；
函数内 _LAST_IDS = ids 只重绑局部名，全局仍是空 → 必须写 _LAST_IDS[:] = ids。

与外部工具对接

本机没有 MAFFT/IQ-TREE 时用纯 Python 跑；有则应交叉验证：

mafft --auto --anysymbol results/msa_trimmed.fasta > results/mafft.aln.fasta
trimal -in results/mafft.aln.fasta -gt 0.9 -fasta > results/mafft.trim.fasta
iqtree2 -s results/mafft.trim.fasta -m LG+G -B 1000 -alrt 1000 -pre results/iqtree -T AUTO

两者拓扑一致 → 结论可靠。不一致时优先信 ML（IQ-TREE），并说明分歧位置。

扩展工作流 A · 基因结构图（外显子/内含子 + 蛋白结构域）

基因家族分析的标准配图之一。当用户说「把序列换成基因结构图」「标注外显子/内含子」
「加上蛋白结构域」时用这条。

数据从哪来（不要用 Ensembl REST，它只覆盖脊椎动物）

水稻基因 ID 要先解析出来（按基因名搜 NCBI Gene 或 Gramene）。Gramene/Ensembl-Plants
的 REST 端点经常不稳定（返回空数组或 503），别把链路建在上面。

解析 Entrez XML 的两个硬要求

1. 标签配对切片，不要用非贪婪正则 —— <Gene-commentary> 是嵌套的，
   正则会跨块（见 pitfalls.md P1-5）。
2. 蛋白长度 = CDS 碱基数 / 3 − 1（CDS 含终止密码子，不减 1 会多算一个氨基酸）。

出图三要素（缺一层读者就得翻别的图）

1. 外显子（盒）/ 内含子（V 形折线，标准记法）/ UTR（矮盒，与 CDS 区分）
2. 蛋白轨道与基因模型同比例绘制 —— 才能同时看到「CDS 相同」与「外显子数不同」
3. 转录方向箭头 —— 负链基因不画箭头，读者必然误读 5'UTR 位置

护栏（绘图前断言，任一不过就中止）

• 外显子数 == NCBI 官方 exon_count property
• 结构域片段 end ≤ protein_length（两个 JSON 的蛋白长度必须一致）
• 外显子不重叠、都在 gene_span 内；外显子长度和 == mRNA 长度
• CDS % 3 == 0；CDS ⊆ 外显子；CDS 不跨内含子

交叉验证（强烈建议做）

拿外显子数去对已发表文献。本项目 OsCATA=4、OsCATC=6 与油棕 CAT 家族论文
及 MDPI 2024 完全吻合 —— 这是对整条取数链路的独立背书，写进报告很有说服力。

扩展工作流 B · 外群定根（判断「谁先分化」）

当用户说「加外群定根」「判断哪个先分化」「无根树要定根」时用这条。

先讲清原理，再动手

• NJ/UPGMA 出的是无根树，无根树上「谁先分化」没有定义 —— 任何边当根都能
  得出不同答案。定根之后才有唯一答案。
• 定根是人为指定根的位置，不是「算出来的」。图注里必须写明。

三个必检项（缺一个结论就不可信）

① 外群够不够远

ratio_min = min(外群→内群距离) / max(内群内部距离)
# ratio_min > 1 才算外群足够远

② 外群是否单系 —— 求外群 LCA，检查其下是否只有外群叶子。非单系 → 定根前提有瑕疵。

③ 多组合稳健性投票 —— 换 4 种外群组合 × 2 种方法，看「谁先分出」是否一致。
只给一棵树是不够的，根的位置依赖外群选择。

分化顺序只用拓扑判断（方法论红线）

「离外群的平均距离」不能用来判断分化先后 —— 距离 = 分化时间 × 进化速率，
两者混淆。正确做法：目标基因组成单系群后，内部最先分出的那一支 = 最早分化。
把这个辨析写进报告和图注，否则会被审稿人抓。见 pitfalls.md P1-4。

外群选择红线

• 必须用单功能过氧化氢酶（clade 1，约 480–515 aa）
• 不能用 katG 类双功能过氧化氢酶-过氧化物酶（约 730 aa，折叠不同）
• 加一个断言：任何序列 > 900 aa → 中止，疑似假基因/融合蛋白

外群非单系时怎么办（重要）

不要为了让图好看而删外群或调整数据。 正确做法是把根因查清并写进论文：

1. 逐条列出每个外群序列的最近邻（在全部序列里）。只要有 1 条的最近邻落在内群，
   外群就不可能单系 —— 这是最直接的证据。
2. 检查外群内部分化深度 vs 外群-内群分化深度。
   若外群内部分化浅于外群-内群分化 → 说明这批外群是谱系特异性旁系同源
   （常来自全基因组复制 WGD），不是「分化前的整体」。
3. 出两幅诊断图（距离矩阵热图 + 最近邻条形图）把根因可视化。
4. 论证结论仍可用的理由：拓扑投票一致性 + 不依赖分支长度 + 与独立证据自洽。

按生物学分组着色

外群和内群要用明显不同的配色，且在图例里按类群分开列（苔藓/藻/真菌/细菌…），
便于审稿人一眼确认外群位置。

交付清单

• [ ] scripts/*.py 可执行脚本（含护栏 + 自检）
• [ ] data/*.fasta 序列（含来源 accession）
• [ ] data/*.json 结构坐标（基因结构图类任务）
• [ ] results/*.fasta 比对（full + trimmed，供外部工具复用）
• [ ] results/distance_matrix_*.tsv
• [ ] results/analysis_report.json 结构化结果
• [ ] figures/*.pdf + *.png 树图 / 基因结构图 / 诊断图
• [ ] README.md 说明参数理由、结论、命名注意事项
• [ ] 方法学文档（分析思路 / 所需数据 / 实施方式 / 结果 / 不确定性）
• [ ] requirements.txt
• [ ] 从零演练：把脚本复制到空目录跑一遍，确认无隐性依赖

> 从零演练是必须做的，不是可选项。 本项目在从零演练时暴露了两个致命 bug
> （--shortcut-dir 参数漏掉、__version__ 改动没落盘），而本地目录因为早有
> 旧文件完全看不出来。演练方式：cp -r scripts data <空目录>，
> 删掉生成的 JSON，按 README 逐步跑通。