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

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

---

## 0. 结论

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

| | 跑前写死的预测 | 实测 | 判定 |
|---|---|---|---|
| **H1** | 洞穴分支上，候选基因判为松弛的比例高于长度配对背景；Fisher p<0.05 且 OR>1 | 候选 **15/84**（17.9%）对背景 **10/84**（11.9%），**OR = 1.61（95% CI 0.68–3.70）**，**Fisher 单侧 p = 0.193**（双侧 0.386） | **被推翻** |
| **H2** | 洞穴分支上，候选基因 K 的中位数低于背景；Mann-Whitney U 双侧 p<0.05 | 中位数 **0.534 对 0.665**（方向对），**U = 3054，双侧 p = 0.133**（单侧 0.067），rank-biserial = 0.13 | **被推翻** |
| **H3**（阴性对照） | 地表分支上不应出现同等富集 | 候选 **6/84**（7.1%）对背景 **8/84**（9.5%），**OR = 0.73（95% CI 0.26–2.17）**，**Fisher 单侧 p = 0.798**；K 中位数 0.775 对 0.912，Mann-Whitney 双侧 p = 0.198 | **成立**（未出现富集，方向还相反） |

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 公开接口，匿名可直接下载。

| 物种 | 登录号 | 注释 | 角色 |
|---|---|---|---|
| *Triplophysa rosa* | GCF_024868665.1 | RefSeq | 洞穴种，前景分支 |
| *Triplophysa dalaica* | GCF_015846415.1 | RefSeq | 地表种，阴性对照的前景分支 |
| *Triplophysa tibetana* | GCA_008369825.1 | 投稿方注释 | 地表种，reference |
| *Barbatula barbatula* | GCF_947034865.1 | RefSeq | 条鳅科外群 |
| *Misgurnus anguillicaudatus* | GCF_027580225.2 | RefSeq | 花鳅科外群 |
| *Danio rerio* | GCF_052040795.1 | RefSeq | 远外群，定根 |

物种树在任何 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 的检出功效随比对长度上升，
不配对时"候选集显著基因更多"完全可以纯由长度差异产生。

| | 数 |
|---|---|
| 计划成对基因 | 136 对 |
| 实际尝试 | 86 对 |
| **进入统计的成对基因** | **84 对（168 个基因）** |
| RELAX 成功拟合 | 341 次 |
| 拟合失败 | 3 次（`pop4` 两条分支、`ccz1` 洞穴分支；整对剔除） |
| CPU 时间 | 115.1 分钟 |

---

## 2. 结果

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

| | 松弛 | 未判为松弛 | 比例 |
|---|---:|---:|---:|
| 候选基因 | **15** | 69 | 17.9% |
| 长度配对背景 | **10** | 74 | 11.9% |

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（洞穴） | p（洞穴） | 地表分支同一基因 |
|---|---|---:|---:|---|
| **cry5** | 隐花色素 / 光受体 | 0.000 | 1.6e-04 | ns |
| **ednrbb** | 色素细胞发育 | 0.132 | 2.8e-04 | **p=0.046，也松弛** |
| oca2 | 黑色素合成 | 0.000 | 2.9e-03 | ns |
| mipb | 晶状体内在蛋白 | 0.194 | 4.8e-03 | ns |
| pmelb | 黑色素体基质 | 0.101 | 5.7e-03 | ns |
| vsx2 | 视网膜发育 | 0.000 | 8.0e-03 | ns |
| pde6c | 视锥光转导 | 0.265 | 9.7e-03 | ns |
| slc24a5 | 色素沉着 | 0.201 | 1.1e-02 | ns |
| six6a | 眼发育 | 0.401 | 1.2e-02 | **p=0.010，也松弛** |
| lim2.2 | 晶状体膜蛋白 | 0.000 | 1.5e-02 | ns |
| gnb3b | 转导蛋白 β | 0.219 | 2.9e-02 | ns |
| rlbp1b | 视黄醛结合蛋白 | 0.199 | 2.9e-02 | ns |
| hlfb | 昼夜节律转录因子 | 0.000 | 3.5e-02 | ns |
| asmt | 褪黑素合成 | 0.222 | 3.9e-02 | ns |
| crybb1 | 晶状体 β-晶体蛋白 | 0.051 | 4.4e-02 | ns |

判为强化（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 中位数（洞穴） | K 中位数（地表） | 配对差中位数 | Wilcoxon 双侧 p |
|---|---:|---:|---:|---:|
| 候选 | 0.534 | 0.775 | −0.189 | 0.114 |
| **背景** | 0.665 | 0.912 | −0.203 | **0.009** |

**洞穴分支上 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` 从跑前落盘起一个字节没动。