204 lines
10 KiB
Markdown
204 lines
10 KiB
Markdown
# 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 会自动回退到：
|
||
1) tests/output/Results/Toxins/Shotter 子目录；若仍不可写，
|
||
2) 当前工作目录下的 ./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) 菌株×目标“目”分数热力图
|
||
```python
|
||
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) 单个菌株：命中×目标“目”贡献热图
|
||
```python
|
||
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` 两列。
|
||
|
||
## 复现命令
|
||
```bash
|
||
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）与图像解读
|
||
|
||
- 生成“菌株×目标目”热力图：
|
||
```bash
|
||
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。
|
||
|
||
- 生成指定菌株的“命中×目标目贡献”热力图：
|
||
```bash
|
||
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）需成对出现，否则降权。
|
||
|
||
1) 单命中相似性权重 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。
|
||
|
||
2) 单命中对各“目”的贡献 c(order)：
|
||
|
||
c(order) = w_hit × P(order | name/subfamily/family)。
|
||
|
||
回退顺序：name → subfamily → family。仍无证据时：
|
||
- 若能解析到家族但无证据 → 贡献计入 other；
|
||
- 若无法解析家族前缀 → 贡献计入 unknown。
|
||
|
||
3) 菌株层合成（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。可在展示时隐藏这两列或单列占比说明。
|