蛋白序列功能注释（证据收集）
何时使用
- 用户丢来一段氨基酸序列或 CDS，问「能得到哪些生物学信息」「有什么结构功能特征」
- 需要判断一个未知蛋白：是不是激酶 / 有没有跨膜段 / 有没有信号肽 / 属于哪个家族
- 需要为某个基因的蛋白产物写功能注释、准备实验设计（突变哪个残基、截短到哪里）
- 需要同源检索、找直系同源、确认物种归属
不适用：只想知道一条序列的分子量（直接手算更快）；需要系统发育树/比对（本 skill 不做进化树，见文末）。
铁律：先确认输入是核酸还是蛋白
用户说「蛋白序列」时，十有八九粘的是 CDS。 实测案例：一条标注为「蛋白质氨基酸序列」
(标注为 RLCK185) 的输入，实际是 1476 nt 的 CDS。
- 判定方法：
is_nucleotide()—— 非 ACGTUN 字符占比 < 15% 且前 60 nt 含ATGG。 - 绝不能默认相信用户的标注。 翻译后必须做 QC：起始
ATG、终止TAA/TAG/TGA、
内部终止子必须为 0、len % 3 == 0。 - 翻译正确性的最终验证：BLAST 打 Swiss-Prot，若首位命中为同一长度、100% 一致度
→ 该物种的参考序列存在，且你的翻译正确。这一步是整条链路的锚点，务必做。
工作流
python scripts/pseq_evidence.py <seq_file> --outdir <DIR> --name <LABEL> --email you@example.com
# 只想跑本地计算（无网/离线）：加 --skip-web
一条命令跑完全部自动化步骤，产出到 --outdir：
| 文件 | 内容 |
|---|---|
protein.txt / protein.faa |
翻译产物 |
evidence_local.json |
组成、理化、KD 扫描、无序、低复杂度、激酶基序、PTM 候选、疏水矩 |
evidence_structure.json |
逐残基 pLDDT、二级结构串、分区 pLDDT 统计、低置信区段 |
ebi_iprscan.tsv / ebi_blast.tsv |
原始结果（保留原始文件，便于复核） |
uniprot_<ACC>.txt |
UniProt 平文件（审阅级功能/定位/PTM/突变体） |
af_model.pdb |
AlphaFold 模型 |
fig_domains.svg / fig_tracks.svg |
结构域图谱 / 多轨剖面 |
DIGEST.txt |
所有数字汇总一处，写报告时直接取用 |
报告怎么写（分工）
脚本只负责收集证据和出图，叙述性报告由你来写。结构建议按用户要求的维度组织：
- 数据溯源与质控（先说明输入是 CDS 还是蛋白，给翻译 QC 表 + BLAST 身份验证）
- 氨基酸组成与理化性质
- 信号肽与跨膜区
- 结构域与功能位点（给"催化亚域逐位定位"表）
- 二级与三级结构（pLDDT 分区 + 二面角二级结构）
- 同源序列与功能注释（给 BLAST 表 + UniProt 功能/GO/互作）
- 翻译后修饰位点（分清"实验验证"与"序列预测候选"两类）
- 基于序列能得到的研究结论（分「可直接下的结论」和「需验证的假说」）
- 方法、数据来源与可复现性（含局限清单）
必须遵守的诚实性检查清单
这些是踩过的坑，不写进报告会得出错误结论：
- 信号肽/跨膜不能声称来自 SignalP/TMHMM。 现行 InterProScan 版本因授权已移除这两个应用，
跑出来不会有它们的条目。用 Kyte-Doolittle 扫描（窗口 19，阈值 1.6）替代，并明确写成
「间接推断」，建议补跑 SignalP 6.0 / DeepTMHMM。 - PTM 分两级写。 UniProt 里有
MOD_RES的才是实验位点；基序扫描出来的全是候选。
若实验位点不能被任何通用基序预测到，这不是脚本 bug —— 反而要写出来，说明它是
"非共有位点（non-consensus）"，只能靠质谱/突变确定。别让工具假装预测到了。 - 基序匹配必须叠加定位合理性。 例子：
N-X-S/T命中了一个 N-糖基化位点，但该蛋白无信号肽、
无跨膜段、定位胞质侧 → 生物学上不可能被糖基化。纯序列扫描不叠加定位信息就会给假阳性。
同理，N 端 G-x-G-x-x-G 常被误判成 ATP 结合环，必须用 PROSITE PS00107 的区间和实验注解校准。 - 低 pLDDT 区的二级结构不可信。 AlphaFold 在无序区常"编"出螺旋。报告里必须点出
"xx–yy 位归为螺旋，但该区 pLDDT < 50，可能是伪影"，不能让读者当成事实。 - 说明二级结构的方法：本 skill 是主链二面角粗分类，不是 DSSP。DSSP 需要实验结构。
- 催化环
YRDvsHRD：不要断言活性受影响。催化碱是不可替换的 Asp，
His 只参与结构；要下结论必须做突变体实验。能说的和不能说的要分开写。 - 功能结论来自文献/UniProt，不是从序列"算"出来的。 序列只能推出"它属于哪类蛋白"，
推不出"它在抗病里起多大作用"。报告里要把这两类结论分开陈述。 - 交叉验证关键数值。
gravy_own与 Biopython 的gravy应完全一致；
MW_Da_own_residue_sum应与 UniProt 的 MW 一致（差 < 1 Da）。不一致就先查代码，
别急着写报告 —— 这类差异往往意味着静默回退（某段代码根本没跑到）。
已知的接口坑（照抄即可，别重新试错）
| 坑 | 正确做法 |
|---|---|
iprscan5 的 stype |
只能是 p / n。传 protein → HTTP 400 |
ncbiblast 的 stype |
只能是 protein / dna。传 p → HTTP 400 |
| 两个服务的词汇不通用 | 同名参数语义不同，一定要看 .../parameterdetails/<name> |
rest.uniprot.org / NCBI 直连 |
部分环境 curl 报 SSL 退出码 35、HTTP 000。用 curl -sk（脚本内已实现 urllib→curl 回退） |
| EBI 可达性自检 | curl -sk -o /dev/null -w "%{http_code}" https://www.ebi.ac.uk/Tools/services/rest/iprscan5/parameters |
| InterProScan 耗时 | 500 aa 蛋白约 1–3 min；脚本按 20 s 轮询，上限 30 min |
| 序列提交 | 用 --data-urlencode "sequence@file" 或 urlencode，别裸拼多行序列 |
| 邮箱 | 必填且必须是合法格式，xxx@example.org 可能被拒；用一个真实格式的地址 |
二级结构二面角规则（如需解释或调参）
α-螺旋 : -100 <= φ <= -30 且 -70 <= ψ <= -5
β-股 : -170 <= φ <= -70 且 (100 <= ψ <= 180 或 -180 <= ψ <= -150)
松散螺旋: -120 <= φ <= -40 且 -60 <= ψ <= 20 -> 归入 H
其余 -> C
预期产物：激酶双叶折叠约 38% H / 28% E / 33% C 量级。
本 skill 不做的事（需要时说清）
- 系统发育树 / 多序列比对：用 MEGA、IQ-TREE，取整家族（含近缘同源）建树
- 复合物对接：AlphaFold-Multimer、HADDOCK
- 表达谱/共表达：UniProt 的
ExpressionAtlas交叉引用、水稻公开转录组 - DSSP 二级结构：需实验结构，拿到 PDB 后另跑
参考基准（用于回归测试本 skill）
拿水稻 OsRLCK185 的 CDS 跑一遍，应得到：
- 1476 nt → 491 aa；BLAST 首位
Q6I5Q6，100% 一致 / 491 aa / E=0 - MW 53 704 Da、pI 9.46（Bjellqvist）、GRAVY −0.610、净电荷 +15
- KD 窗口 19 最大均值 1.284 < 1.6 → 无跨膜螺旋；N 端
MG+ Cys3/Cys6 双酰化 - InterProScan：PROSITE PS50011 80–356；CDD STKc_IRAK；PANTHER PTHR47985:SF16
- 激酶基序：
GEGGFG@87、VAVK@105、YRD@204、DFG@224、APE@251 - AlphaFold：激酶域 pLDDT ≈ 93.6，N 端 ≈ 49.2、C 端 ≈ 44.7（两端为柔性区）
- 实验 PTM：pS240 / pT241 / pT246（由 OsCERK1 磷酸化）
任一项对不上，先怀疑环境（网络/工具版本）而不是序列。
