高原鳅洞穴种的松弛选择检测：方向对上了，但样本量撑不住

这一轮真跑了：NCBI 六份公开基因组注释 → DIAMOND 双向最佳命中求 1:1 直系同源（15267 个） → MAFFT 密码子比对 → 固定物种树 → HyPhy 2.5.101 RELAX 逐基因拟合 341 次 （84 对基因 × 候选/背景 × 洞穴/地表两个前景分支）。纯 CPU 单机，115 分钟。

0. 结论

一句话：预期的方向在数据里出现了，但强度不够，按跑前写死的判据三条预测里两条被推翻； 而被推翻的那两条，把「这条路值不值得继续投」回答得比一个阳性结果更清楚。

H1 和 H2 的点估计都指向预期方向，但两个 p 值都没到 0.05，OR 的置信区间跨过 1。 按跑前判据记为推翻。更贴近实际的读法是：方向正确、量级中等、样本量撑不住—— 既不是"没有效应"，也不是"有效应"。

必须一起读的三条

(1) 这不是全基因组普查。 直系同源这一层是全基因组的（15267 个，全表随包提供）； 但逐基因 RELAX 跑满全基因组约需 170 小时单机 CPU。本轮在同源全集之上跑了一个 跑前预注册、按比对长度 1:1 配对的子集：136 个候选 + 136 个背景，预算 115 分钟用尽即停， 四臂齐全的成对基因 84 对。执行顺序在任何一次 RELAX 开跑之前由固定种子洗牌写死， 截断发生在与结果无关的随机顺序上，不引入选择偏倚。范围差距在跑前写入 claims.json。

(2) 就算跑满 136 对，H1 也不会显著。 按实测比例外推，136 对时 Fisher 单侧 p ≈ 0.115； 在 OR ≈ 1.6 上要达到 80% 功效，每组约需 441 个基因，而这份候选清单在 1:1 同源集里只有 136 个。 瓶颈不在"把剩下的跑完"。

(3) 洞穴种 n=1。 NCBI 上高原鳅属 17 份组装里只有 3 份带注释，洞穴种只有 Triplophysa rosa（GCF_024868665.1）。这条末端分支上的谱系特异效应（有效群体大小、 世代时间、组装与注释质量）与"洞穴适应"完全混在一起，本轮拆不开。 所以本轮结果只能用于假设生成，不能作为"洞穴适应导致视觉基因松弛"的证据。

1. 数据与设计

全部来自 NCBI Datasets v2 公开接口，匿名可直接下载。

物种树在任何 RELAX 生产运行之前固定（659 个背景基因串联密码子 p-distance，NJ）：

`` (Danio_rerio,(Misgurnus_anguillicaudatus,(Barbatula_barbatula,(Triplophysa_rosa,(Triplophysa_dalaica,Triplophysa_tibetana))))); ``

候选集 = 预注册的 242 个视觉 / 光转导 / 晶状体 / 昼夜节律 / 色素基因符号，136 个在 1:1 同源集中存在。 清单按标准功能通路从斑马鱼命名法逐类枚举，没有可引用的外部出处；若要写进论文， 建议改用 GO 定义（GO:0007602 / GO:0007623 / GO:0048066），但换定义会换掉分母。

背景集 = 从 13753 个有正式基因符号的 1:1 同源基因里用固定种子抽 700 个作池， 再按比对长度与每个候选 1:1 最近邻配对。必须长度配对：RELAX 的检出功效随比对长度上升， 不配对时"候选集显著基因更多"完全可以纯由长度差异产生。

2. 结果

2.1 H1：计数上的富集——方向对，不显著

odds ratio = 1.609（Haldane 修正 1.582，95% CI 0.677–3.698），Fisher 单侧 p = 0.193。 15 与 10 相差 5 个基因；几个基因的浮动就能让 OR 在 1.0 到 2.3 之间变，这就是 CI 跨过 1 的含义。

2.2 H2：K 分布位移——同样方向对，不显著

候选 K 中位数 0.534、背景 0.665；Mann-Whitney 双侧 p = 0.133，rank-biserial 0.134（小效应）。 补充一个非预注册、更贴合配对设计的检验：配对内 Wilcoxon 双侧 p = 0.409，支持更弱。 这个数对结论不利，一并列出。

2.3 洞穴分支上判为松弛的候选基因

完整矩阵在 work/results/relax_matrix.tsv。判为松弛（K<1 且 p<0.05）的 15 个候选：

判为强化（K>1 且 p<0.05）的候选 4 个：pax6a、guca1d、cry3a、crx；背景侧松弛 10 个、强化 3 个。

三条限制必须一起读：

1. 单基因 p 未做多重校正。 洞穴分支 168 个基因上按 BH 控制 FDR<0.05，只有 3 个活下来： cry5、ednrbb（候选）和 actr8（背景）。
2. ednrbb、six6a 在地表分支上同样松弛，不是洞穴特有。扣掉后，过 FDR 且洞穴特有的候选只剩 cry5。
3. 多个 K 落在模型边界（K=0.000 是下界），不是精确点估计。

这张表是假设生成清单，不是结论清单。

3. 阴性对照：这把尺子在量什么

3.1 H3 成立

同一批 84 对基因、同一棵树、同一套流程，前景换成地表种 T. dalaica：候选 7.1% 对背景 9.5%， OR = 0.731，Fisher 单侧 p = 0.798。洞穴分支 OR=1.61、地表分支 OR=0.73，计数侧方向相反。但按预注册要求并列的连续量给出的是另一幅图： 地表分支上候选基因的 K 中位数同样低于背景（0.775 对 0.912，差 −0.137，rank-biserial 0.115）， 与洞穴分支（0.534 对 0.665，差 −0.131，rank-biserial 0.134）方向相同、幅度相当。 所以 H3 只在「计数式富集」这一个口径上成立，不能读成「这把尺子是干净的」； 预注册里说可估的「分支×基因集」交互，本轮没有做正式检验。

3.2 同一批对照也暴露了对结论不利的一面

每个基因在两条分支上的 K 配对比较（非预注册）：

洞穴分支上 K 整体偏低，在背景基因上反而更显著。 洞穴分支 168 个基因里名义 p<0.05 的有 32 个 （零假设期望约 8.4 个）——信号存在，但铺在整条分支上，不集中在视觉基因上。 这正是 n=1 混杂在数据上的样子：有效群体大小、注释质量、或全基因组尺度的适应，本轮都排除不掉。

4. 局限（按对结论的威胁排序）

1. 洞穴种 n=1，与谱系特异效应拆不开；结果只能用于假设生成。
2. 样本量不足以支撑显著性；本轮的"不显著"更可能是功效问题，而非"效应不存在"的证据。
3. 候选集可能有偏，方向未定：在洞穴种里已丢失 / 注释不到的视觉基因进不了 1:1 同源集。 实测 12 个候选在地表种和外群里都有、在洞穴种里没有，反方向只有 7 个——但 12 比 7 本身不显著（双侧精确二项检验 p=0.36），两物种质控后的总基因数也本来就差 2.5%（23097 对 23681），而且本轮没有与背景基因的缺失率作对照、没有在基因组上核查假基因化；注释质量差异也能造成同样的不对称。所以只能说这个方向的偏倚可能存在，不能说它已经把松弛程度压低了。
4. 单基因 p 未做多重校正；过 BH FDR 且洞穴特有的候选只有 cry5。
5. RBH 不等于真正的直系同源，未做共线性验证。
6. 注释来源异质：T. tibetana 用投稿方注释，分支长度系统性偏长。
7. RELAX 只检验选择强度的整体变化，不定位位点，不涉及调控区、表达与表型。

1、2、5、6 四条在跑之前写入了 claims.json。

5. 要把它推到有定论，需要什么

• 加大背景而非跑满候选：功效瓶颈在背景侧。按实测比例外推，136 个候选对 1000 个全基因组背景基因， Fisher 单侧 p ≈ 0.043。这是设计计算，不是承诺——比例本身的不确定性很大。
• 把洞穴种从 1 个变成 2–3 个：对另外 14 份无注释的高原鳅属组装做同源引导基因预测， 让前景变成多条独立洞穴分支。这是本轮唯一能被彻底解决的结构性问题。
• 把 cry5 追到位点级（FEL / MEME / branch-site），判断是真的功能松弛还是比对或注释假象。

附：怎么核对每一个数

``bash python3 work/code/08_analyze.py . # 主分析（scipy） python3 work/code/09_recheck_headline.py . # 独立复算（不依赖 scipy） ``

两条实现互不读对方的输出，给出的 Fisher p、odds ratio、K 中位数、完成对数全部一致。 预注册指纹在 work/results/PREREG.sha256，本轮结束时 claims.json 从跑前落盘起一个字节没动。