lingyuzeng/bttoxin-pipeline

Fork 0

Files

zly c9be0694ab Update shotter math documentation with Logit Prior algorithm

2026-01-18 20:34:03 +08:00

34 KiB

Raw Blame History

Bttoxin_Shotter 数学与算法全解（Typora 版）

目标：把 BtToxin_Digger 的命中结果 + BPPRC Specificity Database 的活性数据，转换成
「给定一个菌株，它对哪些昆虫目 / 物种最可能有杀虫活性？」的概率型评分。

本文用 Typora 兼容的公式语法（行内 $...$ ，行间 $$...$$），把所有用到的数学公式、符号和数据来源都说明清楚，并解释为什么这样设计、有哪些假设和局限。

1. 整体思路

用一句话概括 Shotter 的核心想法：

先问“每条毒素命中长得像谁（相似性）？”再问“这些已知毒素在 BPPRC 里主要杀什么虫（特异性）？”，最后把多条证据组合成菌株层的总体风险（noisy-OR）。

更正式一点的步骤：

从 BPPRC 的 toxicity-data.csv 里，统计出：
- 每个 Bt 毒素（名字 / 亚家族 / 家族）
- 对各个昆虫“目”（order）和“物种”（species）的经验活性分布
  → 数学上记作 $P(\text{order} \mid \text{毒素})$、$P(\text{species} \mid \text{毒素})$。
从 BtToxin_Digger 的 All_Toxins.txt 里，取出每条命中的：
- Identity（序列一致性）
- Coverage（比对覆盖度）
- HMM 标记（有没有 HMM 结构域命中） → 折算成一条 0–1 的相似性权重 $w_i$。
对于菌株里的每条命中 $i$，用它 best hit 的毒素名/家族，到第 1 步的索引里查：
- 这条命中，对每个昆虫目/物种的单条贡献：
```
c_i(\text{order}) = w_i \times P(\text{order} \mid \text{该命中的毒素名字/家族})
```
对同一菌株的所有命中，利用 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}} 是该菌株里的全部命中集合。
得到：
- strain_target_scores.tsv：菌株 × 目标“目”的分数矩阵；
- 如果 BPPRC 有物种信息，还会有 strain_target_species_scores.tsv：菌株 × 目标物种矩阵；
- toxin_support.tsv：逐命中 × “目/物种”的贡献矩阵，可用来画 per_hit_C15.png 这种“哪条毒素负责哪个靶标”的可视化。

2. 记号与数据来源一览

2.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$。
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，用死亡率描述毒力。
输出文件
- 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$：

按 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$。
定义毒力权重：
```
\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 在三个层级上建索引：

精确名字：
$\text{name}_r = \text{toxin}_r$，例如 Cry21Aa3；
亚家族：
例如把 Cry21Aa3 截断成 Cry21A；
家族：
再截成 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$：

计算该 key 在各个目上的总权重和：


Z^{\text{order}}_K = \sum_o S_{\text{order}}[K, o]

如果 $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 的设计是：

对每条命中，优先在 name 层级 查 $P(\cdot \mid \text{name})$；
如果该名字没有任何证据（$Z^{\text{order}}_{\text{name}} = 0$），则退到 subfamily 层级；
再不行退到 family 层级；
如果在 family 层级也没有证据，又有两种情况：
- 可以解析出家族前缀（例如知道它是 Cry999 这种新家族），但 BPPRC 里还没有任何活性记录：
  → 这类证据被记到特殊的 other 列中，不去抬高任何具体的“目”。
- 连家族前缀也解析不出来（名字不符合 Bt 命名规则）：
  → 记到 unknown 列中。

这就是 strain_target_scores.tsv 里 other / unknown 的来源。

4. 从 Digger 命中计算相似性权重 `w_i`

4.1 命中是什么？

在 BtToxin_Digger 中，给定一个菌株基因组，流程大致是：

ORF/CDS 预测；
CDS 翻译成蛋白；
用 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 条命中：

从 Digger 的 Hit_ID 得到名字 $\text{name}_i$，解析出对应的：
- $\text{subfamily}_i$，如 Cry21A；
- $\text{family}_i$，如 Cry21。
按以下顺序尝试从特异性索引中取分布：
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 会解析错误。

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 实现逻辑的文字化说明，如需与最新代码完全严格对齐，请以代码仓库中的实现为准。

11. 加入 3 项“基因组背景特征”到菌株×靶标分数（ZWA/Thu/TAA BGC / Mobilome / CRISPR）

你新增的三类指标，本质上都是“对同一个菌株的全基因组背景做先验修正（strain-level priors）”。

Shotter v1 的 \text{Score}(\text{strain},\text{order}) 仅由“毒素命中→靶标分布→noisy-OR”构成（第 3–6 节）。

现在我们引入三类菌株层特征，用来调整（boost / penalize）每个靶标分数，同时保持：

分数仍在 $[0,1]$；

解释仍然清晰：毒素命中给“证据”，背景特征给“先验”。

三类特征都来自序列比对/检测的结果：

ZWA/Thu/TAA：对已知 BGC/关键酶蛋白做 BLAST/HMM（或用 BGC 预测器）→ 得到“存在/不存在”。

Mobilome：对转座酶、质粒复制蛋白、噬菌体标志蛋白等做 BLAST/HMM/注释统计 → 得到数量。

CRISPR/Cas：对 Cas 蛋白做 HMM/BLAST + 阵列检测（repeat/spacer）→ 判定“完整/不完整/不存在”。

11.1 为什么这三项可以影响“杀虫活性潜力”？（三条关键因果链）

你要向专家说明的不是“它们直接提高毒力”，而是：它们改变了菌株在进化与生态上获得/维持杀虫模块的概率。

(A) ZWA/Thu/TAA 三类“杀虫素 BGC”：存在即为直接功能模块证据（正向）

关键因果链（可直接放到讨论里）：

这些基因簇（BGC）代表一种可表达的生物合成路径（酶系/修饰/转运/自抗等模块齐备）。
若用序列比对/结构域检测确认该路径的关键酶与骨架基因存在，则“产物存在的可行性”显著上升。
因此，即使 Cry/Vip 等蛋白毒素命中较少，BGC 的存在也代表一种独立（正交）的杀虫潜力来源。

结论：存在（1）应当加分；不存在（0）不加分。

(B) 移动元件（转座酶/质粒/噬菌体）数量：mobilome 越丰富，越容易获得/重排毒素模块（总体正向，但要饱和）

关键因果链：

Bt 及 B. cereus group 中，大量杀虫相关基因（包括 cry 等）常与质粒/移动 DNA 库强相关，且移动元件（IS/转座子等）参与基因重排、模块拼装与在质粒上的迁移。
mobilome 越丰富，意味着可重排/可迁移的 DNA 元件越多，越容易出现：
- 新毒素模块的获得（HGT/共轭/转导等）；
- 现有毒素模块的复制、重排、组合与剂量效应；
- 数据库未覆盖的新型“other/unknown”杀虫因子的潜力。
但 mobilome 指标也会受组装质量、注释阈值影响，因此应做“边际递减/饱和”，避免被噪声拉爆。

结论：数量越多总体加分（正向），但用饱和函数。

(C) CRISPR/Cas 完整度：越完整，越像“限制外源 DNA 的屏障”，从而降低获得毒素/质粒库的先验（负向）

关键因果链：

CRISPR-Cas 的核心生态功能之一是抵御外源遗传元件（质粒/噬菌体等），本质上会对水平基因转移（HGT）形成选择压力。
在 Bt 中，杀虫谱与毒素库的快速扩展常与质粒/移动 DNA 库相关；
因此，当 CRISPR/Cas 更完整且更可能功能健全时，菌株对外源质粒/移动模块的“进入与稳定维持”通常更困难，导致“获取/更新毒素库”的先验下降。
反过来，CRISPR 缺失/失活的菌株更可能处于“更开放的 mobilome 交换状态”，更容易累积可迁移毒素模块。

结论：对“杀虫潜力总评分”应采用 不存在 > 不完整 > 完整 的单调顺序（即 CRISPR 越完整越压分）。

注：CRISPR 的实际生态效应可能随环境与 anti-CRISPR 等因素复杂化，但作为“全基因组先验”，上述方向能提供最稳定、最可解释的单调修正。

11.2 三类指标的输入格式（对应你的量化规则）

你希望的量化方式：

三种杀虫素（分别为 ZWA、Thu、TAA）生物合成基因簇：存在（1）、不存在（0）。
移动元件（转座酶、质粒、噬菌体）：用数量表示。
CRISPR/Cas 系统：完整（2）、不完整（1）、不存在（0）。

记作：

$b_Z, b_T, b_A \in {0,1}$（对应 ZWA/Thu/TAA）
$m \in \mathbb{N}_0$（mobilome 总计数，或分项加和）
$c \in {0,1,2}$（CRISPR 状态：0=不存在，1=不完整，2=完整）

11.3 为什么不用把这三类特征直接塞进 noisy-OR？（建模选择）

noisy-OR 适合合成“多条毒素命中对同一靶标的独立贡献”。而这三类新特征是“菌株整体背景”，并不对应某个具体命中的 $c_i(\text{order})$。

因此最稳妥的做法是：

先按 v1 算出 $S_\text{tox}(\text{strain},\text{order})$（即第 6 节的 Score）；
再用一个“先验修正项”把它调整为 $S_\text{final}$。

11.4 推荐的并入方式：logit 先验修正（保持 `[0,1]` 且可解释）

令：


S_{\text{tox}}(\text{strain},o) = 1 - \prod_{i\in H_{\text{strain}}}\left[1-c_i(o)\right]

我们定义最终分数：


S_{\text{final}}(\text{strain},o)=\sigma\Big(\operatorname{logit}(S_{\text{tox}}(\text{strain},o)+\varepsilon) + \Delta(\text{strain})\Big)

其中：

\sigma(x)=\dfrac{1}{1+e^{-x}} 是 sigmoid；
$\operatorname{logit}(p)=\ln\dfrac{p}{1-p}$；
\varepsilon 是极小正数（如 $10^{-6}$），避免 \logit(0) 或 $\logit(1)$。

把三项特征写进 $\Delta(\text{strain})$：


\Delta(\text{strain}) =
\beta_Z\,b_Z + \beta_T\,b_T + \beta_A\,b_A
+ \beta_M\,g(m)
+ \beta_C\,h(c)

其中 g(m) 与 h(c) 用来保证“饱和”和“单调方向”。

11.4.1 mobilome 饱和函数 $g(m)$（正向，边际递减）

推荐二选一：

对数饱和：


 g(m)=\ln(1+m)

线性截断：


 g(m)=\min\left(1,\frac{m}{K}\right)

K 取一个经验阈值（例如 50 或 100），取决于你统计的元件类型与注释粒度。

11.4.2 CRISPR 单调映射 $h(c)$（“完整压分”）

为了实现“不存在 > 不完整 > 完整”的顺序，直接把 c\in\{0,1,2\} 映射到 ${1,0.5,0}$：


 h(c)=1-\frac{c}{2}

解释：

不存在 $c=0\Rightarrow h=1$（最大加成）
不完整 $c=1\Rightarrow h=0.5$（中等加成）
完整 $c=2\Rightarrow h=0$（不加成，相当于“最压分”）

只要取 $\beta_C>0$，就得到你想要的确定性方向：

CRISPR 越完整 → h(c) 越小 → \Delta 越小 → S_{\text{final}} 越低

11.5 这三项对评分是“增加还是降低”？（给你一个可写进算法的确定性结论）

在上述形式下，只要你取参数满足：

\beta_Z,\beta_T,\beta_A > 0
\beta_M > 0
\beta_C > 0

则得到确定性结论：

ZWA/Thu/TAA：存在（1）一定增加分数；不存在（0）不增加。
移动元件数量：数量越多一定增加分数（经 g(m) 饱和）。
CRISPR：不存在（0）加成最大；不完整（1）次之；完整（2）最小（最压分）。

11.6 将“三项特征”并入 Shotter 的输出矩阵（菌株×靶标）

你最终仍输出一个矩阵：

strain_target_scores.tsv（或新文件名）里，每个单元格从 S_{\text{tox}} 替换为 $S_{\text{final}}$。

也就是说，Shotter 原本的 per-hit 解释仍然保留（因为 S_{\text{tox}} 没变）；三项新特征只是在菌株层做一个整体“上移/下移”的可解释修正。

11.7 Mermaid：在原流程上加入“背景特征先验”分支（可选）

如果你想把这三项清晰地画进流程图，可以追加一个“Genome context”分支：

flowchart TB
    G["Genome .fna"] --> ORF[ORF/CDS prediction]
    ORF --> TR[Translate CDS to proteins]

    TR --> ID["Toxin identification\n(BLAST/HMM)"]
    ID --> ALL["All_Toxins.txt"]

    BPPRC["BPPRC specificity CSV\n(toxicity-data.csv)"] --> IDX["Build specificity index\nP(order|·), P(species|·)"]

    ALL --> PARSE["Per-hit parsing"]
    IDX --> WHIT["Compute w_i"]
    WHIT --> PHIT["Per-hit scoring\nc_i(order)"]
    PHIT --> COMB["Per-strain noisy-OR\nS_tox(strain,order)"]

    subgraph CONTEXT["Genome context features (BLAST/HMM/detectors)"]
        BGC["ZWA/Thu/TAA BGC\n(0/1)"]
        MGE["Mobilome count m"]
        CRISPR["CRISPR state c\n(0/1/2)"]
    end

    TR --> BGC
    TR --> MGE
    TR --> CRISPR

    COMB --> ADJ["Logit prior adjust\nS_final = sigmoid(logit(S_tox)+Δ)"]
    CONTEXT --> ADJ

    ADJ --> OUT2["strain_target_scores.tsv\n(strain × order, updated)"]

12. 参考文献（用于支撑第 11 节的因果链，建议在论文/汇报中引用）

下面列出的是“因果链”的常用支撑点：

Bt 毒素与质粒/移动 DNA 的关联；

CRISPR-Cas 对质粒/HGT 的屏障作用；

Bt 的移动基因库与致病/宿主谱的关系。

Fiedoruk K. et al. Genetic Environment of cry1 Genes Indicates Their Common Origin. (2017).（cry 基因位于 Bt 质粒并与特定复制系统/遗传环境相关）
Lechuga A. et al. Completed Genomic Sequence of Bacillus thuringiensis… (2020).（Bt 的昆虫致病能力与质粒携带 cry 基因相关的概述）
Thomas D.J.I. et al. Transfer of plasmid pBC16 between Bacillus thuringiensis… (2002).（Bt 大质粒可经类似共轭机制传播，获得 cry 质粒可扩展生态位/宿主）
Gillis A. et al. Role of plasmid plasticity and mobile genetic elements in… (2018).（综述：Bt 移动基因库与致病/生活史相关）
Marraffini L.A. & Sontheimer E.J. CRISPR Interference Limits Horizontal Gene Transfer in Staphylococci. Science (2008).（经典实验证据：CRISPR 干预能限制共轭/转化等 HGT）
Wheatley R.M. & MacLean R.C. CRISPR-Cas systems restrict horizontal gene transfer… (2020).（综述：CRISPR 对质粒/ICE 等 HGT 的普遍限制作用）

34 KiB Raw Blame History Unescape Escape