add math description
This commit is contained in:
685
docs/shotter_math_full_zh_typora.md
Normal file
685
docs/shotter_math_full_zh_typora.md
Normal file
@@ -0,0 +1,685 @@
|
||||
|
||||
# Bttoxin_Shotter 数学与算法全解(Typora 版)
|
||||
|
||||
> 目标:把 **BtToxin_Digger** 的命中结果 + **BPPRC Specificity Database** 的活性数据,转换成
|
||||
> 「给定一个菌株,它对哪些昆虫**目 / 物种**最可能有杀虫活性?」的概率型评分。
|
||||
|
||||
本文用 Typora 兼容的公式语法(行内 `$...$`,行间 `$$...$$`),把所有用到的数学公式、符号和数据来源都说明清楚,并解释为什么这样设计、有哪些假设和局限。
|
||||
|
||||
---
|
||||
|
||||
## 1. 整体思路
|
||||
|
||||
用一句话概括 Shotter 的核心想法:
|
||||
|
||||
> **先问“每条毒素命中长得像谁(相似性)?”再问“这些已知毒素在 BPPRC 里主要杀什么虫(特异性)?”,最后把多条证据组合成菌株层的总体风险(noisy-OR)。**
|
||||
|
||||
更正式一点的步骤:
|
||||
|
||||
1. 从 BPPRC 的 `toxicity-data.csv` 里,统计出:
|
||||
- 每个 Bt 毒素(名字 / 亚家族 / 家族)
|
||||
- 对各个昆虫“目”(order)和“物种”(species)的**经验活性分布**
|
||||
→ 数学上记作 $P(\text{order} \mid \text{毒素})$、$P(\text{species} \mid \text{毒素})$。
|
||||
|
||||
2. 从 BtToxin_Digger 的 `All_Toxins.txt` 里,取出每条命中的:
|
||||
- Identity(序列一致性)
|
||||
- Coverage(比对覆盖度)
|
||||
- HMM 标记(有没有 HMM 结构域命中)
|
||||
→ 折算成一条 0–1 的**相似性权重** $w_i$。
|
||||
|
||||
3. 对于菌株里的每条命中 $i$,用它 best hit 的毒素名/家族,到第 1 步的索引里查:
|
||||
- 这条命中,对每个昆虫目/物种的**单条贡献**:
|
||||
|
||||
$$
|
||||
c_i(\text{order}) = w_i \times P(\text{order} \mid \text{该命中的毒素名字/家族})
|
||||
$$
|
||||
|
||||
4. 对同一菌株的所有命中,利用 **noisy-OR** 公式,在菌株层合成:
|
||||
|
||||
$$
|
||||
\text{Score}(\text{strain}, \text{order}) = 1 - \prod_{i \in H_{\text{strain}}} \left[1 - c_i(\text{order})\right]
|
||||
$$
|
||||
|
||||
其中 $H_{\text{strain}}$ 是该菌株里的全部命中集合。
|
||||
|
||||
5. 得到:
|
||||
- `strain_target_scores.tsv`:菌株 × 目标“目”的分数矩阵;
|
||||
- 如果 BPPRC 有物种信息,还会有 `strain_target_species_scores.tsv`:菌株 × 目标物种矩阵;
|
||||
- `toxin_support.tsv`:逐命中 × “目/物种”的贡献矩阵,可用来画 `per_hit_C15.png` 这种“哪条毒素负责哪个靶标”的可视化。
|
||||
|
||||
---
|
||||
|
||||
## 2. 记号与数据来源一览
|
||||
|
||||
### 2.1 数据文件
|
||||
|
||||
1. **Digger 输出(序列相似性)**
|
||||
- `All_Toxins.txt`(必须):
|
||||
每行是一条命中(hit),包含:
|
||||
- `Strain`:菌株名
|
||||
- `Protein_ID`:该菌株中蛋白的 ID
|
||||
- `Hit_ID`:best hit 的已知 Bt 毒素名字(如 `Cry21Aa3`)
|
||||
- `Identity`:与该毒素的氨基酸一致性(百分数)
|
||||
- `Align_len` / `Query_len` / `Subject_len`:比对长度与两边长度
|
||||
- `HMM` 标记:是否有 HMM 域命中
|
||||
|
||||
在脚本中,这些字段被用来计算 $id_i$、$cov_i$、$hmm_i$。
|
||||
|
||||
2. **BPPRC Specificity Database(经验活性)**
|
||||
- `toxicity-data.csv`(只用阳性记录):
|
||||
- `protein` / `toxin`:实验中使用的 Bt 毒素名,如 `Cry21Aa3`;
|
||||
- `target_order`:目标昆虫目,如 `Rhabditida`、`Lepidoptera` 等;
|
||||
- `target_species`(可选):目标物种学名;
|
||||
- `activity`:是否有活性(只保留“有活性”的记录);
|
||||
- `lc50` + `lc50_unit`:若有 LC50,单位如 `ppm`、`µg/g` 等;
|
||||
- `percentage_mortality` 或文字描述:若无 LC50,用死亡率描述毒力。
|
||||
|
||||
3. **输出文件**
|
||||
- `toxin_support.tsv`:每条命中 $i$ 对每个“目/物种”的贡献 $c_i(\cdot)$;
|
||||
- `strain_target_scores.tsv`:菌株 × 目标“目”的得分 $\text{Score}(\text{strain}, \text{order})$;
|
||||
- `strain_target_species_scores.tsv`:菌株 × 目标物种的得分 $\text{Score}(\text{strain}, \text{species})$;
|
||||
- `strain_scores.json`:上述结果的 JSON 版本。
|
||||
|
||||
### 2.2 数学符号
|
||||
|
||||
为了后面讲公式方便,统一一下符号:
|
||||
|
||||
- $\text{strain}$:一个菌株,例如 `C15`;
|
||||
- $H_{\text{strain}}$:这个菌株中所有命中的集合;
|
||||
- $i \in H_{\text{strain}}$:第 $i$ 条命中(hit);
|
||||
- $\text{order}$:某个昆虫目,如 `Lepidoptera`、`Diptera` 等;
|
||||
- $\text{species}$:某个具体物种(如果 CSV 有);
|
||||
- $\text{name}$:某个毒素的精确名字,例如 `Cry21Aa3`;
|
||||
- $\text{subfamily}$:亚家族记号,例如 `Cry21A`;
|
||||
- $\text{family}$:家族记号,例如 `Cry21`。
|
||||
|
||||
对应到 Digger / BPPRC 数据:
|
||||
|
||||
- 对第 $i$ 条命中:
|
||||
- $\text{name}_i$:Digger 里 `Hit_ID` 的值;
|
||||
- $\text{family}_i$:从 `Hit_ID` 解析出来的家族前缀(例如把 `Cry21Aa3` 截成 `Cry21`);
|
||||
- $id_i$:Identity / 100,落在 $[0,1]$;
|
||||
- $cov_i$:比对覆盖度,典型定义是
|
||||
|
||||
$$
|
||||
cov_i = \frac{\text{Align\_len}_i}{\text{参考毒素长度}}
|
||||
$$
|
||||
|
||||
- $hmm_i$:若 HMM 有命中,取 $1$,否则 $0$。
|
||||
|
||||
- 对 BPPRC 的第 $r$ 行实验记录:
|
||||
- $\text{toxin}_r$:毒素名;
|
||||
- $\text{order}_r$、$\text{species}_r$:该条记录对应的目标目/物种;
|
||||
- $\text{LC50}_r$、单位 $\text{unit}_r$;
|
||||
- $\text{mortality}_r$:死亡率或文字,如 “80%”、“some”、“stunting”。
|
||||
|
||||
---
|
||||
|
||||
## 3. 从 BPPRC 构建特异性索引 $P(\text{order} \mid \cdot)$
|
||||
|
||||
### 3.1 只用阳性记录
|
||||
|
||||
首先,在 `toxicity-data.csv` 中只保留 `activity` 表示“有活性”(Yes / positive 等)的记录。理由:
|
||||
|
||||
- 阴性记录经常语义不清晰(剂量太低?实验条件不合适?);
|
||||
- v1 版本采用保守策略:**只用“反复证明打得动”的正证据**,不去尝试建一个完整的剂量-反应模型。
|
||||
|
||||
所以从现在开始,所有行 $r$ 都是 **阳性活性记录**。
|
||||
|
||||
### 3.2 每条记录的“毒力权重” $\text{potency}_r$
|
||||
|
||||
我们希望不同文献、不同剂量单位下的结果,能压缩成一个 0–1 的权重:$\text{potency}_r \in [0,1]$。
|
||||
直觉是:“LC50 越低 / 死亡率越高 → 毒力越强 → 权重越高”。
|
||||
|
||||
#### 3.2.1 按单位桶分组
|
||||
|
||||
因为 LC50 单位多种多样,不能直接混在一起比较:
|
||||
|
||||
- 将记录按 `lc50_unit` 分成若干“单位桶”,例如:
|
||||
- 饮食掺喂实验中 `ppm` ≈ `µg/g` ≈ `ug/g` ≈ `μg/g`,合并到同一桶;
|
||||
- 其它单位(如 `µg/cm^2`, `mg/L` 等)各自成桶。
|
||||
|
||||
在同一个桶里,单位一致,才能比较大小。
|
||||
|
||||
#### 3.2.2 有 LC50 的记录:反分位归一
|
||||
|
||||
对于某个单位桶 $u$,设其中有 $n_u$ 条记录,集合记作 $R_u$。
|
||||
对每条记录 $r \in R_u$:
|
||||
|
||||
1. 按 LC50 从小到大排序,给它一个“从 0 到 1 的排名百分位”:
|
||||
|
||||
$$
|
||||
\text{rank\_pct}_r =
|
||||
\begin{cases}
|
||||
0, & n_u = 1 \\
|
||||
\dfrac{\text{rank}_u(r)-1}{n_u - 1}, & n_u > 1
|
||||
\end{cases}
|
||||
$$
|
||||
|
||||
- 这里 $\text{rank}_u(r)$ 是 LC50 排序后的位次(最小的 LC50 rank = 1)。
|
||||
- 因此 LC50 最小的记录 $\text{rank\_pct}=0$,最大的为 $1$。
|
||||
|
||||
2. 定义毒力权重:
|
||||
|
||||
$$
|
||||
\text{potency}_r = 1 - \text{rank\_pct}_r
|
||||
$$
|
||||
|
||||
- LC50 最低(毒力最强) ⇒ $\text{potency}=1$;
|
||||
- LC50 最高(毒力最弱) ⇒ $\text{potency}=0$;
|
||||
- 中间线性过渡。
|
||||
|
||||
**为什么不用直接做 $\text{potency} = 1 / \text{LC50}$?**
|
||||
|
||||
- 不同单位之间本来就不能直接比;
|
||||
- 同一单位桶内部的相对排序比绝对数值更可靠;
|
||||
- 分位数归一对极端值不敏感,适合文献数据这种“噪声大 + 范围广”的情况。
|
||||
|
||||
#### 3.2.3 没有 LC50,只有死亡率的记录
|
||||
|
||||
如果某条记录没有 LC50,但有死亡率或文字描述,则按照经验规则赋值:
|
||||
|
||||
- 死亡率 ≥ 80% 或文字为 “greater than 80%” →
|
||||
|
||||
$$
|
||||
\text{potency}_r = 0.90
|
||||
$$
|
||||
|
||||
- 50–80%(包括 60–80%、60–100% 等区间)→
|
||||
|
||||
$$
|
||||
\text{potency}_r = 0.65
|
||||
$$
|
||||
|
||||
- 0–60% / 0–50% / “25%” / “some” / “stunting”等弱效应 →
|
||||
|
||||
$$
|
||||
\text{potency}_r = 0.25
|
||||
$$
|
||||
|
||||
这些常数(0.90 / 0.65 / 0.25)不是“统计拟合”的结果,而是**人工选择的分级:强 / 中 / 弱**。
|
||||
|
||||
#### 3.2.4 既没有 LC50 也没有死亡率文本
|
||||
|
||||
如果一条记录只告诉你“有活性”,但没有剂量信息:
|
||||
|
||||
$$
|
||||
\text{potency}_r = 0.55
|
||||
$$
|
||||
|
||||
解释:
|
||||
|
||||
- 介于“中等毒力”(0.65)和“最低档有毒”(0.25)之间;
|
||||
- 既体现“确实有活性”的正证据,又不过分抬高其权重。
|
||||
|
||||
### 3.3 从单条记录到毒素级别的总权重
|
||||
|
||||
每条记录 $r$ 还包含一个毒素名 $\text{toxin}_r$、目标目 $\text{order}_r$ 和可选的目标物种 $\text{species}_r$。
|
||||
|
||||
Shotter 在三个层级上建索引:
|
||||
|
||||
1. **精确名字**:
|
||||
$\text{name}_r = \text{toxin}_r$,例如 `Cry21Aa3`;
|
||||
2. **亚家族**:
|
||||
例如把 `Cry21Aa3` 截断成 `Cry21A`;
|
||||
3. **家族**:
|
||||
再截成 `Cry21`。
|
||||
|
||||
对于每一个层级的 key(记为 $K$,可以是 name / subfamily / family),分别计算:
|
||||
|
||||
- 在“目标目”维度的总毒力:
|
||||
|
||||
$$
|
||||
S_{\text{order}}[K, o] =
|
||||
\sum_{\substack{r \text{ 属于阳性记录} \\
|
||||
\text{key}(r) = K \\
|
||||
\text{order}_r = o}}
|
||||
\text{potency}_r
|
||||
$$
|
||||
|
||||
- 如果 CSV 有物种信息,“目标物种”维度的总毒力:
|
||||
|
||||
$$
|
||||
S_{\text{species}}[K, s] =
|
||||
\sum_{\substack{r \text{ 属于阳性记录} \\
|
||||
\text{key}(r) = K \\
|
||||
\text{species}_r = s}}
|
||||
\text{potency}_r
|
||||
$$
|
||||
|
||||
直觉:
|
||||
|
||||
> 每一行实验是一个“该毒素对某个目标有怎样强度的正证据”,
|
||||
> 把同一毒素、同一目标的所有实验权重加起来,就得到**该毒素对该目标的总支持度**。
|
||||
|
||||
### 3.4 归一化成概率分布 $P(\text{order} \mid K)$
|
||||
|
||||
对每一个 key $K$:
|
||||
|
||||
1. 计算该 key 在各个目上的总权重和:
|
||||
|
||||
$$
|
||||
Z^{\text{order}}_K = \sum_o S_{\text{order}}[K, o]
|
||||
$$
|
||||
|
||||
2. 如果 $Z^{\text{order}}_K > 0$,则定义:
|
||||
|
||||
$$
|
||||
P(\text{order}=o \mid K) =
|
||||
\frac{S_{\text{order}}[K, o]}{Z^{\text{order}}_K}
|
||||
$$
|
||||
|
||||
对所有 order 求和等于 1。
|
||||
|
||||
同理,如果有物种信息:
|
||||
|
||||
$$
|
||||
P(\text{species}=s \mid K) =
|
||||
\frac{S_{\text{species}}[K, s]}
|
||||
{\sum_{s'} S_{\text{species}}[K, s']}
|
||||
$$
|
||||
|
||||
这一步的含义是:
|
||||
|
||||
> “**在所有已报道的靶标中,这个毒素的毒力是如何在不同‘目/物种’之间分布的?**”
|
||||
|
||||
例如,对某个毒素名字 $K = \text{Cry21Aa3}$,可能得到类似分布:
|
||||
|
||||
- $P(\text{Rhabditida} \mid K) = 0.8$
|
||||
- $P(\text{Thysanoptera} \mid K) = 0.2$
|
||||
- 其它目为 0。
|
||||
|
||||
### 3.5 回退(back-off)与 `other` / `unknown`
|
||||
|
||||
现实中很多毒素只有家族层面数据,比如只知道 `Cry21`,没有具体到 `Cry21Aa3`;
|
||||
因此 Shotter 的设计是:
|
||||
|
||||
1. 对每条命中,优先在 **name 层级** 查 $P(\cdot \mid \text{name})$;
|
||||
2. 如果该名字没有任何证据($Z^{\text{order}}_{\text{name}} = 0$),则退到 **subfamily 层级**;
|
||||
3. 再不行退到 **family 层级**;
|
||||
4. 如果在 family 层级也没有证据,又有两种情况:
|
||||
- 可以解析出家族前缀(例如知道它是 `Cry999` 这种新家族),但 BPPRC 里还没有任何活性记录:
|
||||
→ 这类证据被记到特殊的 `other` 列中,不去抬高任何具体的“目”。
|
||||
- 连家族前缀也解析不出来(名字不符合 Bt 命名规则):
|
||||
→ 记到 `unknown` 列中。
|
||||
|
||||
这就是 `strain_target_scores.tsv` 里 `other` / `unknown` 的来源。
|
||||
|
||||
---
|
||||
|
||||
## 4. 从 Digger 命中计算相似性权重 $w_i$
|
||||
|
||||
### 4.1 命中是什么?
|
||||
|
||||
在 BtToxin_Digger 中,给定一个菌株基因组,流程大致是:
|
||||
|
||||
1. ORF/CDS 预测;
|
||||
2. CDS 翻译成蛋白;
|
||||
3. 用 BLAST + HMM 去已知 Bt 毒素数据库里做毒素识别。
|
||||
|
||||
`All_Toxins.txt` 中的每一行就是一条 **命中(hit)**:
|
||||
|
||||
- 这条蛋白 “像” 某个已知 Bt 毒素($\text{name}_i$);
|
||||
- Digger 给出 Identity、比对长度、HMM 命中等信息;
|
||||
- Shotter 需要把“像的程度”压缩成 0–1 权重 $w_i$。
|
||||
|
||||
### 4.2 Identity、Coverage 与 HMM 变量
|
||||
|
||||
对第 $i$ 条命中:
|
||||
|
||||
- $id_i$:Digger 的 `Identity` 字段除以 100,变成 0–1 小数;
|
||||
- $cov_i$:比对覆盖度,典型取值为
|
||||
|
||||
$$
|
||||
cov_i = \frac{\text{Align\_len}_i}{\text{参考毒素长度}_i}
|
||||
$$
|
||||
|
||||
- $hmm_i$:如果 HMM 检测到相应家族的结构域,则为 1,否则为 0。
|
||||
|
||||
### 4.3 base 分段函数:家族 / 亚家族阈值
|
||||
|
||||
Bt 毒素命名规范中,**家族 / 亚家族** 的划分大致用 Identity 0.45 和 0.78 作为分界。
|
||||
Shotter 按此定义一个分段函数:
|
||||
|
||||
$$
|
||||
\text{base}_i =
|
||||
\begin{cases}
|
||||
1, & id_i \ge 0.78 \text{ 且 } cov_i \ge 0.80 \\
|
||||
\dfrac{id_i - 0.45}{0.78 - 0.45}, & 0.45 \le id_i < 0.78 \\
|
||||
0, & \text{否则}
|
||||
\end{cases}
|
||||
$$
|
||||
|
||||
含义:
|
||||
|
||||
- 当 $id_i \ge 0.78$ 且 $cov_i \ge 0.80$ 时:
|
||||
- 认为这条命中和参考毒素至少属于**同一个亚家族**,高度可靠,给 $\text{base}_i = 1$;
|
||||
- 当 $0.45 \le id_i < 0.78$:
|
||||
- 认为像是同家族但亚家族不确定,于是 base 在 0–1 之间线性上升;
|
||||
- 当 $id_i < 0.45$:
|
||||
- 序列过于遥远,直接视为无效命中,$\text{base}_i = 0$。
|
||||
|
||||
coverage 条件 $cov_i \ge 0.80$ 是为了过滤掉“只覆盖一小截结构域”的伪阳性。
|
||||
|
||||
### 4.4 合并 coverage 与 HMM 加成,得到 $w_i$
|
||||
|
||||
最终的相似性权重定义为:
|
||||
|
||||
$$
|
||||
w_i = \min\left(1,\ \text{base}_i \cdot cov_i + 0.1 \cdot hmm_i\right)
|
||||
$$
|
||||
|
||||
解释:
|
||||
|
||||
- $\text{base}_i \cdot cov_i$:
|
||||
base 体现 Identity 的可信度,再乘以 coverage,避免只有一小段高 identity 时权重过高。
|
||||
- $0.1 \cdot hmm_i$:
|
||||
如果 HMM 有命中($hmm_i=1$),则额外 $+0.1$,表示“结构域层面的再确认”;
|
||||
但上限由 $\min(1,\dots)$ 限制在 1 以内。
|
||||
|
||||
**为什么 HMM 只加 0.1?**
|
||||
|
||||
- HMM 命中说明有“家族结构域”存在,是一种很好的辅助证据;
|
||||
- 但 HMM 也可能比较宽泛,不必高于完整序列比对;
|
||||
- 因此把它当作“尾数加分”,而不是主导因子。
|
||||
|
||||
### 4.5 配对家族(partner families)的下调
|
||||
|
||||
某些 Bt 毒素家族只有在“配对”时才有完整毒力,例如:
|
||||
|
||||
- Vip1 / Vip2
|
||||
- Vpa / Vpb
|
||||
|
||||
Shotter 的处理是:
|
||||
|
||||
- 对于在“需要配对”的家族集合中的一条命中 $i$:
|
||||
- 检查同一菌株中是否存在“配对家族”的命中 $j$,且 $w_j \ge 0.3$;
|
||||
- 如果找不到这样的搭档,则**下调当前命中的权重**:
|
||||
|
||||
$$
|
||||
w_i \leftarrow 0.2 \times w_i
|
||||
$$
|
||||
|
||||
- 下调而不直接置零,是因为:
|
||||
- 仍然保留“弱证据”,比如可能另一条 partner 被漏检;
|
||||
- 但在没有配对证据的情况下,不应让这种命中主导靶标预测。
|
||||
|
||||
---
|
||||
|
||||
## 5. 单命中对各目标目/物种的贡献 $c_i(\cdot)$
|
||||
|
||||
有了相似性权重 $w_i$ 和特异性分布 $P(\text{order} \mid K)$ / $P(\text{species} \mid K)$,每条命中就可以被看成是对不同靶标的“分布式投票”。
|
||||
|
||||
### 5.1 为每条命中选择合适的 key(name / subfamily / family)
|
||||
|
||||
对第 $i$ 条命中:
|
||||
|
||||
1. 从 Digger 的 `Hit_ID` 得到名字 $\text{name}_i$,解析出对应的:
|
||||
- $\text{subfamily}_i$,如 `Cry21A`;
|
||||
- $\text{family}_i$,如 `Cry21`。
|
||||
2. 按以下顺序尝试从特异性索引中取分布:
|
||||
1. 若存在 $P(\text{order} \mid \text{name}_i)$,则取 $K_i = \text{name}_i$;
|
||||
2. 否则若存在 $P(\text{order} \mid \text{subfamily}_i)$,则 $K_i = \text{subfamily}_i$;
|
||||
3. 否则若存在 $P(\text{order} \mid \text{family}_i)$,则 $K_i = \text{family}_i$;
|
||||
4. 若三者都不存在:
|
||||
- 若能解析出 family,则认为“这是一个有家族但 BPPRC 没记录的新毒素”,贡献放入 `other`;
|
||||
- 若连 family 也解析不出,则放入 `unknown`。
|
||||
|
||||
### 5.2 贡献公式:$c_i(\text{order}) = w_i \times P(\text{order} \mid K_i)$
|
||||
|
||||
对每个命中 $i$ 和每个目标目 $\text{order}$,定义:
|
||||
|
||||
$$
|
||||
c_i(\text{order}) = w_i \times P(\text{order} \mid K_i)
|
||||
$$
|
||||
|
||||
解释:
|
||||
|
||||
- $P(\text{order} \mid K_i)$:经验上,**这个毒素的毒力主要集中在哪些目**;
|
||||
- $w_i$:这条命中到底“有多像”这个毒素;
|
||||
- 乘积 $w_i \times P(\text{order} \mid K_i)$ 就是:
|
||||
|
||||
> “在考虑序列相似度之后,这条命中对某个目有多强的杀虫活性支持度”。
|
||||
|
||||
对 `other` / `unknown` 的处理:
|
||||
|
||||
- 如果 $K_i$ 在任何层级都没有 BPPRC 证据,则:
|
||||
- 若能解析家族前缀:
|
||||
- $c_i(\text{other}) = w_i$,其它目为 0;
|
||||
- 若连家族前缀也没有:
|
||||
- $c_i(\text{unknown}) = w_i$,其它目为 0。
|
||||
|
||||
### 5.3 物种维度:$c_i(\text{species})$
|
||||
|
||||
如果 `toxicity-data.csv` 里有 `target_species` 信息,则同样构建 $P(\text{species} \mid K)$,然后:
|
||||
|
||||
$$
|
||||
c_i(\text{species}) = w_i \times P(\text{species} \mid K_i)
|
||||
$$
|
||||
|
||||
这一部分会被写入 `toxin_support.tsv` 的物种列,并最终合成 `strain_target_species_scores.tsv`。
|
||||
|
||||
---
|
||||
|
||||
## 6. 菌株层合成:noisy-OR 公式
|
||||
|
||||
单条命中的贡献 $c_i(\cdot)$ 只能告诉我们:
|
||||
|
||||
> “如果只有这一条毒素存在,那么它对某个目有多大概率有活性?”
|
||||
|
||||
但实际菌株中往往有多条毒素命中,需要把它们合起来。
|
||||
|
||||
### 6.1 数学形式
|
||||
|
||||
对某个菌株 `strain`,它的命中集合为 $H_{\text{strain}}$。
|
||||
对任意一个目标目 $\text{order}$,定义:
|
||||
|
||||
$$
|
||||
\text{Score}(\text{strain}, \text{order}) =
|
||||
1 - \prod_{i \in H_{\text{strain}}} \left[1 - c_i(\text{order})\right]
|
||||
$$
|
||||
|
||||
同理,如果看物种维度:
|
||||
|
||||
$$
|
||||
\text{Score}(\text{strain}, \text{species}) =
|
||||
1 - \prod_{i \in H_{\text{strain}}} \left[1 - c_i(\text{species})\right]
|
||||
$$
|
||||
|
||||
### 6.2 公式中每个符号的含义
|
||||
|
||||
- $c_i(\text{order})$:第 $i$ 条命中对该目有活性的概率 / 支持度,来源于上一节;
|
||||
- $1 - c_i(\text{order})$:对应地,是“**第 $i$ 条命中对该目不起作用**”的概率;
|
||||
- $\prod_{i \in H_{\text{strain}}} [\dots]$:大写 Π 表示**连乘号**,
|
||||
$\prod_{i} x_i$ 意思是 $x_1 \times x_2 \times \cdots$;
|
||||
- 所以
|
||||
|
||||
$$
|
||||
\prod_{i \in H_{\text{strain}}} [1 - c_i(\text{order})]
|
||||
$$
|
||||
|
||||
可以理解为:
|
||||
|
||||
> “菌株里所有毒素都对这个目不起作用”的概率(假设独立)。
|
||||
|
||||
- 最外面的 $1 - (\cdot)$ 就是:
|
||||
|
||||
> “至少有一条毒素对这个目起作用”的概率。
|
||||
|
||||
这就是经典的 **noisy-OR** 组合方式:
|
||||
逻辑上“OR”(只要至少有一个 cause 生效),但每个 cause 都带有“噪声”(概率 < 1)。
|
||||
|
||||
### 6.3 直观小例子
|
||||
|
||||
假设某个菌株里,对 `Rhabditida` 有两条命中:
|
||||
|
||||
- 命中 1:$c_1 = 0.4$
|
||||
- 命中 2:$c_2 = 0.6$
|
||||
|
||||
那么:
|
||||
|
||||
$$
|
||||
\text{Score} =
|
||||
1 - (1 - 0.4) \times (1 - 0.6)
|
||||
= 1 - 0.6 \times 0.4
|
||||
= 1 - 0.24
|
||||
= 0.76
|
||||
$$
|
||||
|
||||
可以看到:
|
||||
|
||||
- $\text{Score}(\text{strain}, \text{Rhabditida}) = 0.76$ 比单条 0.6 / 0.4 都大;
|
||||
- 但仍然在 0–1 之间,不会像简单相加那样变成 1.0 以上。
|
||||
|
||||
如果有 3 条贡献分别为 0.3、0.3、0.3:
|
||||
|
||||
$$
|
||||
\text{Score} = 1 - 0.7^3 \approx 1 - 0.343 = 0.657
|
||||
$$
|
||||
|
||||
说明多个中等证据叠加起来,也能给出一个比较高的综合得分。
|
||||
|
||||
### 6.4 为什么认为“独立”是合理近似?
|
||||
|
||||
noisy-OR 的核心假设是:不同毒素对同一靶标的“失效”是近似独立的。
|
||||
|
||||
现实中毒素之间可能互相协同或者拮抗,但当前阶段:
|
||||
|
||||
- 没有足够数据建完整的“多毒素交互模型”;
|
||||
- 又希望有一个简单、可解释、单调的组合方式。
|
||||
|
||||
在这样的权衡下,**noisy-OR 是一个非常常见、保守且易于解释的选择**:
|
||||
|
||||
- 任意增加新的命中,只会让 Score 不减(增加或持平),符合常识;
|
||||
- 多条弱证据可以叠加成中等甚至强证据;
|
||||
- 强证据不会被简单叠加到 >1,而是趋近于 1。
|
||||
|
||||
---
|
||||
|
||||
## 7. 从分数到“最大概率靶标”
|
||||
|
||||
有了 `strain_target_scores.tsv` 之后,对任意菌株 `strain`:
|
||||
|
||||
- 可以找到:
|
||||
|
||||
$$
|
||||
\text{best\_order}(\text{strain}) =
|
||||
\arg\max_{\text{order}} \text{Score}(\text{strain}, \text{order})
|
||||
$$
|
||||
|
||||
作为**单一“最大概率”目标目**;
|
||||
|
||||
- 也可以设置阈值,如:
|
||||
- Score ≥ 0.7:强候选靶标;
|
||||
- 0.4–0.7:中等可信度;
|
||||
- < 0.4:弱证据。
|
||||
|
||||
物种层面的 `strain_target_species_scores.tsv` 也可以同样处理,得到:
|
||||
|
||||
- 最可能的靶标物种;
|
||||
- 以及“该菌株在某个目内部,更偏向于打哪些代表物种”。
|
||||
|
||||
---
|
||||
|
||||
## 8. 常见质疑与回答
|
||||
|
||||
这一节列出一些“审稿人/同行可能会问”的问题和解释,方便在 PPT 或论文中使用。
|
||||
|
||||
### Q1:为什么只用阳性记录,不用阴性记录?
|
||||
|
||||
- 阴性结果的语义模糊:
|
||||
“没看到活性”可能是剂量太低、实验条件不合适、靶标抗性等多种原因;
|
||||
- 数据中大量阴性条目缺乏剂量信息,很难和阳性条目放在一个统一框架里比较;
|
||||
- v1 设计选择保守策略:**只累积“肯定打得动”的正证据**。
|
||||
|
||||
### Q2:LC50 / 死亡率是不同单位和实验条件,为什么可以直接比较?
|
||||
|
||||
- LC50 只在**同一单位桶内**被排序和归一化,不跨单位;
|
||||
- 采用**反分位数**而不是直接用 $1/\text{LC50}$,是为了:
|
||||
- 用“谁比谁强”而不是“绝对值差多少”;
|
||||
- 对极端值不敏感;
|
||||
- 不同实验条件确实会带来噪声,但在大量文献叠加时,**稳定的“毒素→靶标”关系**仍然会在分布中凸显出来。
|
||||
|
||||
### Q3:0.45 / 0.78 这两个 Identity 阈值怎么来的?
|
||||
|
||||
- 它们来自 Bt 毒素命名的经验规则:
|
||||
通行做法大致是:
|
||||
- ≥0.95:同一毒素;
|
||||
- ~0.78:亚家族边界;
|
||||
- ~0.45:家族边界。
|
||||
- Shotter 把 0.45–0.78 区间视为“同家族但亚家族不确定”,用线性插值连接 0 和 1;
|
||||
- 这是与 Bt 命名体系兼容的**连续近似**:
|
||||
- Identity 明显高于 0.78 → 强证据;
|
||||
- 明显低于 0.45 → 弱证据,基本忽略;
|
||||
- 中间区域在两个极端之间平滑过渡。
|
||||
|
||||
### Q4:为什么用 noisy-OR 而不是简单相加 / 取最大值?
|
||||
|
||||
- 简单相加:
|
||||
- 多条中等证据累加可能 >1,不再有概率意义;
|
||||
- 不能很好表达“至少有一条毒素起作用”的逻辑。
|
||||
- 取最大值:
|
||||
- 会完全忽略其它毒素的贡献;
|
||||
- 无法体现“多条中等证据也值得重视”。
|
||||
- noisy-OR 的优点:
|
||||
- 单调:增加证据只会增加或保持分数;
|
||||
- 有明确概率解释:$1 - \prod(1-c_i)$;
|
||||
- 实现简单,可读性好,方便向非数学背景使用者解释。
|
||||
|
||||
### Q5:BPPRC 数据不全面,新家族或少量报道会不会被低估?
|
||||
|
||||
- 对于**完全没有 BPPRC 记录的家族**,命中会被计入 `other` 或 `unknown`,不会错误地抬高某个具体目;
|
||||
- 这意味着:
|
||||
- 新家族 / 数据稀少家族 在 v1 中只能作为“待挖掘候选”,而不是高置信度靶标;
|
||||
- 可以在未来版本中补充其它来源(专利、未发表数据、实验室结果)来充实特异性索引。
|
||||
|
||||
---
|
||||
|
||||
## 9. Mermaid 流程图(Typora 友好版)
|
||||
|
||||
> 注意:含有点号或 `<br/>` 的节点需要加引号,否则 Mermaid 会解析错误。
|
||||
|
||||
```mermaid
|
||||
flowchart TB
|
||||
G["Genome .fna"] --> ORF[ORF/CDS prediction]
|
||||
ORF --> TR[Translate CDS to proteins]
|
||||
TR --> ID["Toxin identification<br/>(BLAST/HMM)"]
|
||||
|
||||
subgraph Digger["BtToxin_Digger (Docker)"]
|
||||
ORF
|
||||
TR
|
||||
ID
|
||||
end
|
||||
|
||||
ID --> LIST["*.list / *.gbk"]
|
||||
ID --> ALL["All_Toxins.txt"]
|
||||
|
||||
BPPRC["BPPRC specificity CSV<br/>(toxicity-data.csv)"] --> IDX["Build specificity index<br/>P(order|·), P(species|·)"]
|
||||
|
||||
LIST --> PARSE["Per-hit parsing<br/>(one row per toxin hit)"]
|
||||
ALL --> PARSE
|
||||
|
||||
PARSE --> WHIT["Compute similarity weight w_i"]
|
||||
IDX --> WHIT
|
||||
|
||||
WHIT --> PHIT["Per-hit scoring<br/>c_i(order), c_i(species)"]
|
||||
PHIT --> COMB["Per-strain combination<br/>noisy-OR over hits"]
|
||||
|
||||
COMB --> OUT1["toxin_support.tsv<br/>(hit × order/species contributions)"]
|
||||
COMB --> OUT2["strain_target_scores.tsv<br/>(strain × order scores)"]
|
||||
COMB --> OUT3["strain_target_species_scores.tsv<br/>(strain × species scores, optional)"]
|
||||
COMB --> OUT4["strain_scores.json"]
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## 10. 使用建议
|
||||
|
||||
- 做 **方法学 PPT / 论文方法部分** 时:
|
||||
- 可直接引用第 3–6 节的公式,配合 1–2 个数值小例子;
|
||||
- Mermaid 图可以直接贴到 Typora、GitHub 或 Obsidian 中渲染。
|
||||
- 给 **非专业听众** 讲解时:
|
||||
- 重点讲 “像谁 + 那个家族历史上主要打什么虫 + 多条证据用 noisy-OR 合成” 这条主线;
|
||||
- 正文公式保留 score 这一条即可。
|
||||
|
||||
_本文是对 Bttoxin_Shotter v1 实现逻辑的文字化说明,如需与最新代码完全严格对齐,请以代码仓库中的实现为准。_
|
||||
Reference in New Issue
Block a user