lingyuzeng/bttoxin-pipeline

Fork 0

Files

hotwa 82cfb13e50 add usage markdown file

2025-11-21 20:23:15 +08:00

10 KiB

Raw Blame History

Bttoxin_Shoter 结果文件说明与可视化建议

本页介绍 Shotter 的输出文件结构、字段含义、解读建议，以及如何快速画热力图等可视化。

文件一览（默认写入 tests/output/Results/Toxins 或可写回退目录）

toxin_support.tsv
- 逐命中（per-hit）对各目标“目”的贡献矩阵。
strain_target_scores.tsv
- 每个菌株（strain）对各“目”的综合分数矩阵（0–1）。
strain_scores.json
- 与上表相同信息的 JSON 版。

注：All_Toxins.txt 是 BtToxin_Digger 的输出，Shotter 的输入；toxin_support.tsv/strain_target_scores.tsv/strain_scores.json 是 Shotter 的新输出文件。

若默认目录不可写（如容器生成目录归 root 所有），Shotter 会自动回退到：

tests/output/Results/Toxins/Shotter 子目录；若仍不可写，
当前工作目录下的 ./shotter_outputs。命令行也可用 --output_dir 指定任意可写目录。

字段与含义

toxin_support.tsv（per-hit 矩阵）

Strain：菌株名（来自 All_Toxins.txt 的第一列）
Protein_id：Digger 报告的蛋白/ORF 标识
Hit_id：Digger 命中名（去掉 “-other” 尾巴后的标准名）
Identity：序列一致性（0–1），由 Digger 的 Identity 除以 100 后得到
Aln_length / Hit_length：比对长度与参考长度
Coverage：覆盖度（Aln_length/Hit_length，0–1）
HMM：Digger 的 HMM 列（YES/NO）
Family：解析到的家族键（如 Cry1、Vip3）；无法解析则为 unknown
NameKey：用于回退检索的名称键（通常等于 Hit_id）
PartnerFulfilled：是否满足配对家族（如 Vip1/Vip2、Vpa/Vpb）
Weight：该命中的相似性权重 w_hit，公式：
- 若 id≥0.78 且 cov≥0.80：base=1.0；elif 0.45≤id<0.78：base=(id-0.45)/(0.78-0.45)；else base=0
- w_hit = min(1, base×coverage + 0.1×HMM)
- 若需配对但未满足：w_hit×=0.2
后续每一列：一个 target_order（如 Lepidoptera、Coleoptera… 以及 other、unknown），值为该命中对该“目”的贡献
- 贡献 c(order)=w_hit×P(order|name/subfamily/family)
- 若能解析出家族但 CSV 无证据 → 计入 other
- 若无法解析家族 → 计入 unknown

strain_target_scores.tsv（per-strain 矩阵）

Strain：菌株名
每个 target_order 一列（含 other、unknown）：
- Score(strain, order)=1−∏(1−c(order))（noisy-OR 合成）
- 落在 0–1。值越大，综合证据越强。可自定阈值作“高置信”标记（例如 ≥0.5）。

strain_scores.json

与 strain_target_scores.tsv 等价的 JSON，便于程序消费。

如何解读

per-hit：适合审查具体哪条毒蛋白命中贡献了哪些“目”。Weight 与 Coverage/Identity/HMM 协同体现可信度；PartnerFulfilled=NO 时表示该家族可能需要配对，单独命中仅做弱证据。
per-strain：纵览一个菌株对哪些“目”潜在活性更强。可以取 Top-N order 查看最可能靶标，也可设置阈值筛“高置信”目标。
列 other/unknown：
- other：可解析到家族，但 CSV 中无对应靶标证据（回退失败）
- unknown：无法从命中名中解析出家族前缀

热力图可视化示例

1) 菌株×目标“目”分数热力图

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

df = pd.read_csv("tests/output/Results/Toxins/strain_target_scores.tsv", sep="\t")
# 如果程序回退到 ./shotter_outputs，请改成相应路径

orders = [c for c in df.columns if c != "Strain"]
mat = df.set_index("Strain")[orders]

plt.figure(figsize=(0.6*len(orders)+3, 0.4*len(mat)+2))
sns.heatmap(mat, cmap="viridis", vmin=0, vmax=1, cbar_kws={"label": "score"})
plt.title("Strain vs Target Orders (Shotter)")
plt.tight_layout()
plt.show()

2) 单个菌株：命中×目标“目”贡献热图

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

hits = pd.read_csv("tests/output/Results/Toxins/toxin_support.tsv", sep="\t")
# 选择一个菌株
s = "97-27"
sub = hits.query("Strain == @s").copy()
orders = [c for c in sub.columns if c not in [
    "Strain","Protein_id","Hit_id","Identity","Aln_length","Hit_length","Coverage",
    "HMM","Family","NameKey","PartnerFulfilled","Weight"
]]
mat = sub.set_index("Hit_id")[orders]

plt.figure(figsize=(0.6*len(orders)+3, 0.4*len(mat)+2))
sns.heatmap(mat, cmap="magma", vmin=0, vmax=1, cbar_kws={"label": "contrib"})
plt.title(f"Per-hit contributions for {s}")
plt.tight_layout()
plt.show()

常见问题

为什么不再做 BLAST？
- Shotter 直接使用 Digger 的相似性结果与 HMM 标记；避免重复计算，提高速度与可复现性。
为什么会有 PermissionError？
- 某些 Digger 输出目录可能由容器创建（root 所有），普通用户不可写。可传入 --output_dir 指向可写目录，或使用程序的自动回退（Shotter 子目录或 ./shotter_outputs）。
是否只看 Top1？
- 矩阵保留了所有“目”的得分，更利于发现“次 high”的潜在靶标。若需要，我可以在输出中追加 TopOrder/TopScore 两列。

复现命令

python scripts/bttoxin_shoter.py \
  --toxicity_csv Data/toxicity-data.csv \
  --all_toxins tests/output/Results/Toxins/All_Toxins.txt \
  --output_dir tests/output/Results/Toxins
# 若目录不可写可改为： --output_dir tests/output/Results/Toxins/Shotter

绘图脚本用法（CLI）与图像解读

生成“菌株×目标目”热力图：
```
python scripts/plot_shotter.py \
  --strain_scores shotter_outputs/strain_target_scores.tsv \
  --out_dir shotter_outputs
```
- 行=菌株（Strain）。
- 列=目标“目”（含 other/unknown）。
- 颜色=综合分数 Score(strain, order) ∈ [0,1]。分数越大，综合证据越强。非“能杀”，而是“潜在活性支持度”。建议论文中给出阈值（如≥0.5）或标注Top-N order。
生成指定菌株的“命中×目标目贡献”热力图：
```
python scripts/plot_shotter.py \
  --strain_scores shotter_outputs/strain_target_scores.tsv \
  --toxin_support shotter_outputs/toxin_support.tsv \
  --out_dir shotter_outputs \
  --per_hit_strain 97-27
```
- 行=命中蛋白（Hit_id；可在图注中标注对应 Family）。
- 列=目标“目”。
- 颜色=贡献 c(order)=w_hit×P(order|·)。可直观看到“哪个命中倾向于哪个‘目’”。

提示：若想在图中暂时不展示 other/unknown，可在绘图前用 pandas 过滤掉相应列。

命中识别与打分算法（适合论文引用）

记：

id ∈ [0,1]：Identity/100；cov ∈ [0,1]：Coverage=Aln_length/Hit_length；hmm∈{0,1}：HMM 标记。
P(order | name/subfamily/family)：由 BPPRC 阳性数据聚合归一得到的“靶标目”分布。
伙伴家族（如 Vip1/Vip2、Vpa/Vpb）需成对出现，否则降权。

单命中相似性权重 w_hit：

若 id ≥ 0.78 且 cov ≥ 0.80：base = 1.0；否则若 0.45 ≤ id < 0.78：base = (id − 0.45) / (0.78 − 0.45)；否则 base = 0。

w_hit = min(1, base × cov + 0.1 × hmm)。

若命中所属家族需配对但菌株内未检测到配对家族（且其权重≥0.3），则 w_hit := 0.2 × w_hit。

单命中对各“目”的贡献 c(order)：

c(order) = w_hit × P(order | name/subfamily/family)。

回退顺序：name → subfamily → family。仍无证据时：

若能解析到家族但无证据 → 贡献计入 other；
若无法解析家族前缀 → 贡献计入 unknown。

菌株层合成（noisy-OR）：

Score(strain, order) = 1 − ∏_hits (1 − c(order))。

理由与阈值说明：0.45/0.78 阈值源于 Bt 命名规范（家族/亚家族）；cov≥0.8 保护部分匹配；HMM 加0.1为结构域级支持；noisy-OR 可在[0,1]内自然合成多条证据、避免重复计数。论文中可将 Score 视作“相对支持度/潜在概率分数”，并配置信心阈值或校准方案。

other/unknown 的来源与处理

other：
- 能解析到家族（如 Cry1），但在特异性索引（BPPRC CSV）中没有该 name/subfamily/family 的靶标证据条目，于是回退失败，贡献计入 other。
- 也可能该命名属于非经典 Bt 杀虫蛋白分支、或尚无公开的活性证据。
unknown：
- 无法从 Hit_id 解析出家族前缀（例如 InhA2-other、Bmp1-other、Zwa6-other 等不在 Cry/Cyt/Vip/Tpp/Mpp…前缀集合中）。此类命中多为其他功能蛋白或未纳入本索引体系的类别。
为什么不能仅用命名系统直接映射靶标？
- 命名主要反映系统发育/结构域类型；同一（亚）家族中不同变体对受体/虫阶的谱系可能显著不同，且跨“目”的活性也常见。以经验规则硬编码会产生系统性偏差，因此我们采用从文献/专利汇总的阳性证据构建 P(order|·) 的统计方法，更稳健、可更新。
绘图实践：若 other/unknown 占比高，表示大量命中缺乏到“目”级证据或无法解析；可在展示时隐藏这两列，或在方法学中单独说明其来源与占比。

两类热力图在科研中的用途建议

菌株×“目”热力图（per-strain）：适合横向比较多个菌株，展示总体“靶标目”潜在谱系；可配 Top-N 或阈值标注。
命中×“目”热力图（per-hit）：适合阐释某一菌株的组成性证据与机制依据，说明哪些命中推动了特定“目”的得分。
论文建议：两者结合使用——正文给出 per-strain 总览、补充材料给出关键样本的 per-hit 解析图。

示例图解读（常见疑问）

非零表示“有一定证据支持潜在活性”，并非确证致死。建议报告高于阈值（如≥0.5）的“高置信目标目”，并在方法中给出阈值与依据。
“per-hit 图纵坐标是什么？”
- 纵轴是命中蛋白的 Hit_id（可在图注中标注对应 Family）。
- 行=命中蛋白（Hit_id；可在图注中标注对应 Family）。
- 列=目标“目”。
- 颜色=贡献 c(order)=w_hit×P(order|·)。可直观看到“哪个命中倾向于哪个‘目’”。
“why other 很多？”
- 多因命中属于未纳入家族前缀集合或 CSV 无证据（例如 InhA2/Bmp1/Zwa6 等），因此贡献汇入 other/unknown。可在展示时隐藏这两列或单列占比说明。

10 KiB Raw Blame History Unescape Escape