基因组数据集#
首先，下载了梁老师给我的数据集，我补充了侧柏的基因组
/HDD3/jch/Gymnosperm_genome_26_5_11/genome_data再补充了白豆杉的基因组，但是白豆杉的格式是gbff，搞明白了，现在大家都知道gbff不好用，会选择上传一份新的数据到figshare上
对于gff和fasta不匹配的巨柏基因组，校准了fasta文件
perl -pe 's/^>\S+(.*OriSeqID=)(\S+)/>$2$1$2/' Cgigantea.fa > Cgigantea_renamed.fa && mv Cgigantea_renamed.fa Cgigantea.fa对于地中海柏木也是一样
perl -pe 'if(/^>/){$i++;$_=">Scaffold_".sprintf("%05d",$i)."\t".substr($_,1)}' Csempervirens.fa > Csempervirens_renamed.fa && mv Csempervirens_renamed.fa Csempervirens.fa这俩物种都是2n=11
检测这些数据集有没有去除最长转录本
#!/bin/bash
# 用法: bash check_transcript_filter.sh species_index.tsv 数据目录
TSV="${1:-species_index.tsv}"
DATADIR="${2:-.}"
echo -e "物种\tGFF基因\tGFF转录本\t转录本/基因\tPEP条数\t判定"
tail -n +2 "$TSV" | while IFS=$'\t' read -r name cat latin fa gff pep cds status; do
# 跳过注释和不完整物种
[[ "$name" =~ ^# ]] && continue
[ "$pep" = "-" ] && continue
[ "$gff" = "-" ] && continue
gff_path="$DATADIR/$gff"
pep_path="$DATADIR/$pep"
[ -f "$gff_path" ] || { echo -e "$name\tGFF不存在\t-\t-\t-\t-"; continue; }
[ -f "$pep_path" ] || { echo -e "$name\tPEP不存在\t-\t-\t-\t-"; continue; }
genes=$(awk '$3=="gene"' "$gff_path" | wc -l)
mrnas=$(awk '$3=="mRNA" || $3=="transcript"' "$gff_path" | wc -l)
pep_count=$(grep -c "^>" "$pep_path")
# 基因数为0时特殊处理
if [ "$genes" -eq 0 ]; then
genes=$(awk -F'[=;]' '{for(i=1;i<=NF;i++) if($i=="gene_id") print $(i+1)}' "$gff_path" | sort -u | wc -l)
fi
# 判定
if [ "$mrnas" -eq 0 ] && [ "$genes" -gt 0 ]; then
if [ "$pep_count" -le $((genes + genes/10)) ]; then
verdict="✓ PEP≈基因"
else
verdict="⚠ 需检查"
fi
elif [ "$mrnas" -gt 0 ]; then
diff_gene=$((pep_count - genes))
diff_mrna=$((pep_count - mrnas))
[ $diff_gene -lt 0 ] && diff_gene=$((-diff_gene))
[ $diff_mrna -lt 0 ] && diff_mrna=$((-diff_mrna))
if [ $diff_gene -le $((genes / 50)) ] || [ $diff_gene -le 50 ]; then
verdict="✓ 已筛选"
elif [ $diff_mrna -le $((mrnas / 50)) ] || [ $diff_mrna -le 50 ]; then
verdict="✗ 未筛选"
else
verdict="⚠ 需手动检查"
fi
else
verdict="⚠ 无gene/mRNA"
fi
[[ "$pep" =~ longest ]] && verdict="$verdict [pre-filtered]"
if [ "$genes" -gt 0 ]; then
ratio=$(awk "BEGIN {printf \"%.2f\", $mrnas/$genes}")
else
ratio="N/A"
fi
echo -e "$name\t$genes\t$mrnas\t$ratio\t$pep_count\t$verdict"
donespecies_index.tsv文件
物种名 分类 完整拉丁名 基因组(fa) GFF PEP CDS 完整性
Aalba [G] Abies alba Aalba.fa Aalba.gff Aalba.pep Aalba.cds 完整
Agramineus [AG?] 待确认 Agramineus.fa Agramineus.gff Agramineus.pep Agramineus.cds 完整
Aspinulosa [OT?] 待确认 (疑似 Alsophila spinulosa 桫椤) Aspinulosa.fa Aspinulosa.gff Aspinulosa.pep Aspinulosa.cds 完整
Athaliana [AG] Arabidopsis thaliana Athaliana.fa Athaliana.gff Athaliana.pep Athaliana.cds 完整
Atrichopoda [AG] Amborella trichopoda Atrichopoda.fa Atrichopoda.gff Atrichopoda.pep Atrichopoda.cds 完整
Cgigantea [G] Cupressus gigantea Cgigantea.fa Cgigantea.gff Cgigantea.pep Cgigantea.cds 完整
Cjaponica [G] Cryptomeria japonica Cjaponica.fa Cjaponica.gff Cjaponica.pep Cjaponica.cds 完整
Cpanzhihuaensis [G] Cycas panzhihuaensis Cpanzhihuaensis.fa Cpanzhihuaensis.gff Cpanzhihuaensis.pep Cpanzhihuaensis.cds 完整
Csempervirens [G] Cupressus sempervirens Csempervirens.fa Csempervirens.gff Csempervirens.pep Csempervirens.cds 完整
Fhodginsii [G] Fokienia hodginsii Fhodginsii.fa Fhodginsii.gff Fhodginsii.pep Fhodginsii.cds 完整
Gbiloba [G] Ginkgo biloba Gbiloba.fa Gbiloba.gff Gbiloba.pep Gbiloba.cds 完整
Gmontanum [G] Gnetum montanum Gmontanum.fa Gmontanum.gff Gmontanum.pep Gmontanum.cds 完整
Lkaempferi [G] Larix kaempferi Lkaempferi.fa Lkaempferi.gff Lkaempferi.longest.pep Lkaempferi.longest.cds 完整
Mglyptostroboides [G] Metasequoia glyptostroboides Mglyptostroboides.fa Mglyptostroboides.gff Mglyptostroboides.pep Mglyptostroboides.cds 完整
Mpolymorpha [OT] Marchantia polymorpha Mpolymorpha.fa Mpolymorpha.gff Mpolymorpha.pep Mpolymorpha.cds 完整
Ncolorata [AG] Nymphaea colorata Ncolorata.fa Ncolorata.gff Ncolorata.pep Ncolorata.cds 完整
Osativa [AG] Oryza sativa Osativa.fa Osativa.gff Osativa.pep Osativa.cds 完整
Pabies [G] Picea abies Pabies.fa Pabies.gff Pabies.pep Pabies.cds 完整
Pdensiflora [G] Pinus densiflora Pdensiflora.fa Pdensiflora.gff Pdensiflora.pep Pdensiflora.cds 完整
Pglauca [G] Picea glauca Pglauca.fa Pglauca.gff Pglauca.pep Pglauca.cds 完整
Plambertiana [G] Pinus lambertiana Plambertiana.fa Plambertiana.gff Plambertiana.pep Plambertiana.cds 完整
Platycladus_orientalis [G] Platycladus orientalis Platycladus_orientalis.fa Platycladus_orientalis_mRNA_longest.gff3 Platycladus_orientalis.pep Platycladus_orientalis.cds完整
Pmenziesii [G] Pseudotsuga menziesii Pmenziesii.fa Pmenziesii.gff Pmenziesii.pep Pmenziesii.cds 完整
Pseudotaxus_chienii [G] Pseudotaxus chienii Pseudotaxus_chienii_pep_chr.fa - - - 仅基因组
Ptabuliformis [G] Pinus tabuliformis Ptabuliformis.fa Ptabuliformis.gff Ptabuliformis.pep Ptabuliformis.cds 完整
Ptaeda [G] Pinus taeda Ptaeda.fa Ptaeda.gff Ptaeda.pep Ptaeda.cds 完整
Sgiganteum [G] Sequoiadendron giganteum Sgiganteum.fa Sgiganteum.gff Sgiganteum.pep Sgiganteum.cds 完整
Ssempervirens [G] Sequoia sempervirens Ssempervirens.fa Ssempervirens.gff Ssempervirens.pep Ssempervirens.cds 完整
Tchinensis [G] Taxus chinensis Tchinensis.fa Tchinensis.gff Tchinensis.pep Tchinensis.cds 完整
Tdistichum [G] Taxodium distichum Tdistichum.fa Tdistichum.gff Tdistichum.pep Tdistichum.cds 完整
Tgrandis [G] Torreya grandis Tgrandis.fa Tgrandis.gff Tgrandis.pep Tgrandis.cds 完整
Twallichiana [G] Taxus wallichiana Twallichiana.1.fa Twallichiana.1.gff Twallichiana_longest.pep Twallichiana_longest.cds 完整
Vvinifera [AG] Vitis vinifera Vvinifera.fa Vvinifera.gff Vvinifera.pep Vvinifera.cds 完整
Wmirabilis [G] Welwitschia mirabilis Wmirabilis.fa Wmirabilis.gff Wmirabilis.pep Wmirabilis.cds 完整结果
(SNP_calling) zhanglab@zhanglab-Precision-3680:/HDD3/jch/Gymnosperm_genome_26_5_11/genome_data$ bash check_transcript_filter.sh
物种 GFF基因 GFF转录本 转录本/基因 PEP条数 判定
Aalba 94205 98227 1.04 97750 ✗ 未筛选
Agramineus 13986 13986 1.00 13986 ✓ 已筛选
Aspinulosa 67831 71488 1.05 67831 ✓ 已筛选
Athaliana 27655 48456 1.75 27654 ✓ 已筛选
Atrichopoda 26846 26846 1.00 26846 ✓ 已筛选
Cgigantea 35384 40910 1.16 40910 ✗ 未筛选
Cjaponica 55246 55246 1.00 55246 ✓ 已筛选
Cpanzhihuaensis 0 32353 N/A 32353 ✓ 未筛选
#攀枝花苏铁是通过了的，他gff里就没有基因这一行
#Csempervirens 670001 720048 1.07 365410 ⚠ 需手动检查
#？？？？？？
#明白了，Csempervirens这个物种的注释有问题，不可用
Fhodginsii 50521 50521 1.00 50521 ✓ 已筛选
Gbiloba 27836 27836 1.00 22152 ⚠ 需手动检查
#银杏怎么还能少了一些转录本的？！莫非手动去除了？有可能
Gmontanum 0 27354 N/A 27354 ✓ 未筛选
#这个没问题，可以使用，是因为统计的gene行没有出现。
Lkaempferi 45826 55814 1.22 42608 ⚠ 需手动检查 [pre-filtered]
Mglyptostroboides 32174 32174 1.00 32174 ✓ 已筛选
Mpolymorpha 19287 24674 1.28 19287 ✓ 已筛选
Ncolorata 24105 44024 1.83 24059 ✓ 已筛选
Osativa 42189 52424 1.24 42189 ✓ 已筛选
Pabies 70736 0 0.00 66632 ✓ PEP≈基因
#为什么pep文件还能比基因组上的注释要少啊
Pdensiflora 44233 4637 0.10 44233 ✓ 已筛选
Pglauca 0 0 N/A 567 ⚠ 无gene/mRNA
#这个物种不使用，注释很烂，组装应该也一般。
Plambertiana 38518 0 0.00 38518 ✓ PEP≈基因
#这个物种不使用，注释很烂，肯定没有达到可用的水平。
Platycladus_orientalis 52897 52897 1.00 52897 ✓ 已筛选
Pmenziesii 51419 0 0.00 51419 ✓ PEP≈基因
Ptabuliformis 80495 144584 1.80 144552 ✗ 未筛选
#100%是把俩单倍型都合并了
Ptaeda 51751 0 0.00 51751 ✓ PEP≈基因
#2017年发表的，不可用
Sgiganteum 41631 0 0.00 41631 ✓ PEP≈基因
#这个可以用，张仁纲用了
Ssempervirens 118906 0 0.00 118906 ✓ PEP≈基因
#这个张仁纲用了，不对。
Tchinensis 44770 44770 1.00 44770 ✓ 已筛选
Tdistichum 42307 42307 1.00 44010 ⚠ 需手动检查
#！！！！！！！！！！！！！！
Tgrandis 47089 47089 1.00 47089 ✓ 已筛选
#要换成李俞鹏ba版本的
Twallichiana 37766 47722 1.26 37766 ✓ 已筛选 [pre-filtered]
Vvinifera 31845 55564 1.74 31845 ✓ 已筛选
Wmirabilis 0 26990 N/A 26990 ✗ 未筛选
#这个是没问题的。拯救油松的十几万个mRNA的办法：
#略AGAT去冗余
# 从 AGAT 处理过的 GFF 提取 mRNA ID
awk '$3=="mRNA"' Ptabuliformis.gff | grep -oP 'ID=\K[^;]+' > agat_ids.txt
# seqkit 按 ID 提取
seqkit grep -f agat_ids.txt Ptabuliformis.pep > Ptabuliformis_longest.pep
# 验证
echo "GFF mRNA 数: $(wc -l < agat_ids.txt)"
echo "提取 PEP 数: $(grep -c '^>' Ptabuliformis_longest.pep)"
(SNP_calling) zhanglab@zhanglab-Precision-3680:/HDD3/jch/Gymnosperm_genome_26_5_11/genome_data$ echo "GFF mRNA 数: $(wc -l < agat_ids.txt)"
echo "提取 PEP 数: $(grep -c '^>' Ptabuliformis_longest.pep)"
GFF mRNA 数: 80495
提取 PEP 数: 80464
提取 CDS 数 80464目前能做到蛋白数接近需求数目(去除最长转录本后的数量)，但是这三十多个基因我确实没招了。
白豆杉的基因组文件重命名
awk '
/^>/ {
# 提取第一个空格后的全部内容作为描述（保留原始格式）
space_pos = index($0, " ")
if (space_pos > 0) {
rest = substr($0, space_pos) # 包含前导空格
} else {
rest = ""
}
# 定位 “ZQBX ” 之后的标记部分（到逗号为止）
marker = ""
if (match($0, /ZQBX /)) {
start_pos = RSTART + RLENGTH
# 在 start_pos 之后查找逗号
after = substr($0, start_pos)
comma_in_after = index(after, ",")
if (comma_in_after == 0) {
marker = after
} else {
marker = substr(after, 1, comma_in_after - 1)
}
}
# 去掉 marker 的前导空格
sub(/^ +/, "", marker)
# 根据模式生成新 ID
if (marker ~ /^chromosome /) {
# 提取数字部分（去掉 “chromosome ”）
num = marker
sub(/^chromosome /, "", num)
new_id = "chr" num
} else if (marker ~ /^Scaffold/) {
new_id = marker
} else {
# 兜底：保留原始 accession（如 CM132175.1）
tmp = substr($0, 2) # 去掉 “>”
sub(/ .*/, "", tmp)
new_id = tmp
}
print ">" new_id rest
next
}
{ print }
' Pseudotaxus_chienii.fa > Pseudotaxus_chienii_renamed.faseqkit seq -w 60 Pseudotaxus_chienii_renamed.fa -o Pseudotaxus_chienii_renamed_formatted.fa最后保留下来的物种
(base) zhanglab@zhanglab-Precision-3680:/HDD3/jch/Gymnosperm_genome_26_5_11/genome_data/genome_saved_26_6_4$ cat 21species_working_set.tsv
物种名 分类 完整拉丁名 中文名 基因组(fa) GFF PEP CDS 完整性 用途
# === 核心裸子植物 (16) ===
Cgigantea [G] Cupressus gigantea 巨柏 Cgigantea.fa Cgigantea.gff Cgigantea.pep Cgigantea.cds 完整 内群
Cjaponica [G] Cryptomeria japonica 日本柳杉 Cjaponica.fa Cjaponica.gff Cjaponica.pep Cjaponica.cds 完整 内群
Cpanzhihuaensis [G] Cycas panzhihuaensis 攀枝花苏铁 Cpanzhihuaensis.fa Cpanzhihuaensis.gff Cpanzhihuaensis.pep Cpanzhihuaensis.cds 完整 内群
Fhodginsii [G] Fokienia hodginsii 福建柏 Fhodginsii.fa Fhodginsii.gff Fhodginsii.pep Fhodginsii.cds 完整 内群
Gbiloba [G] Ginkgo biloba 银杏 Gbiloba.fa Gbiloba.gff Gbiloba.pep Gbiloba.cds 完整 内群
Mglyptostroboides [G] Metasequoia glyptostroboides 水杉 Mglyptostroboides.fa Mglyptostroboides.gff Mglyptostroboides.pep Mglyptostroboides.cds 完整 内群
Pdensiflora [G] Pinus densiflora 赤松 Pdensiflora.fa Pdensiflora.gff Pdensiflora.pep Pdensiflora.cds 完整 内群
Platycladus_orientalis [G] Platycladus orientalis 侧柏 Platycladus_orientalis.fa Platycladus_orientalis_mRNA_longest.gff3 Platycladus_orientalis.pep Platycladus_orientalis.cds 完整 内群
Pseudotaxus_chienii [G] Pseudotaxus chienii 白豆杉 Pseudotaxus_chienii.fa Pseudotaxus_chienii.gff Pseudotaxus_chienii.pep Pseudotaxus_chienii.cds 完整 内群
Ptabuliformis [G] Pinus tabuliformis 油松 Ptabuliformis.fa Ptabuliformis.gff Ptabuliformis.pep Ptabuliformis.cds 完整 内群
Sgiganteum [G] Sequoiadendron giganteum 巨杉 Sgiganteum.fa Sgiganteum.gff Sgiganteum.pep Sgiganteum.cds 完整 内群
Tchinensis [G] Taxus chinensis 红豆杉 Tchinensis.fa Tchinensis.gff Tchinensis.pep Tchinensis.cds 完整 内群
Tdistichum [G] Taxodium distichum 落羽杉 Tdistichum.fa Tdistichum.gff Tdistichum.pep Tdistichum.cds 完整 内群
Tgrandis [G] Torreya grandis 香榧 female_chr_Nc.fna female_ba_last.gff3 Tgrandis.pep Tgrandis.cds 完整 内群
Twallichiana [G] Taxus wallichiana 喜马拉雅红豆杉 Twallichiana.1.fa Twallichiana.1.gff Twallichiana_longest.pep Twallichiana_longest.cds 完整 内群
Clanceolata [G] Cunninghamia lanceolata 杉木 Clanceolata.fa Clanceolata.gff Clanceolata.pep Clanceolata.cds 完整 内群
# === 买麻藤类 (2, 裸子植物外类群) ===
Gmontanum [G] Gnetum montanum 买麻藤 Gmontanum.fa Gmontanum.gff Gmontanum.pep Gmontanum.cds 完整 外类群
Wmirabilis [G] Welwitschia mirabilis 百岁兰 Wmirabilis.fa Wmirabilis.gff Wmirabilis.pep Wmirabilis.cds 完整 外类群
# === 被子植物/苔藓外类群 (3) ===
Atrichopoda [AG] Amborella trichopoda 无油樟 Atrichopoda.fa Atrichopoda.gff Atrichopoda.pep Atrichopoda.cds 完整 外类群
Mpolymorpha [OT] Marchantia polymorpha 地钱 Mpolymorpha.fa Mpolymorpha.gff Mpolymorpha.pep Mpolymorpha.cds 完整 外类群
Ncolorata [AG] Nymphaea colorata 蓝睡莲 Ncolorata.fa Ncolorata.gff Ncolorata.pep Ncolorata.cds 完整 外类群
# === 补充信息 ==
# 这个数据集中的相当部分是梁老师下载的，我补充了其中的杉木，侧柏，白豆杉，对于所有物种进行了最长转录本检查。
# 检查了pep数目是否与注释数目符合。对于不符合的物种进行了校正。Python复写代码
我让AI复写了这个代码之后，他把这个代码分成了五部分，包括以下内容：
主要步骤：
读取控制文件(ctl file)，获取物种信息、BLAST程序路径等配置。
执行BLASTP搜索(如需要)，或跳过使用已有结果。
解析BLASTP输出，识别候选嵌合基因的长拷贝(long copy)和短拷贝(short copy)。
根据系统发育树和物种分组(内群/中群/外群)对候选基因进行过滤和展示。
生成HTML格式的最终结果表格。
#!/usr/bin/env python3
import os
import sys
import re
import argparse
import subprocess
from ete3 import Tree
from Bio import SeqIO
def parse_args_and_build_config():
parser = argparse.ArgumentParser(description='GFFRIN: fusion gene analysis')
# 参数定义保持不变
parser.add_argument('--pep_folder', required=True, help='蛋白序列文件夹路径')
parser.add_argument('--tree', required=True, help='系统发育树文件 (Newick格式)')
parser.add_argument('--ingroups', required=True, help='内群物种，逗号分隔')
parser.add_argument('--outgroups', required=True, help='外群物种，逗号分隔')
parser.add_argument('--output_dir', required=True, help='输出目录')
parser.add_argument('--num_threads', type=int, default=10, help='BLASTP 线程数，默认 10')
args = parser.parse_args()
# 1. 仅保存树文件路径，不做任何解析
tree_path = args.tree
if not os.path.isfile(tree_path):
print(f"树文件不存在: {tree_path}")
sys.exit(1)
# 2. 扫描 PEP 文件夹，构建 protein_path
protein_path = {}
found_species = [] # 替代原来的 intree_check
if not os.path.isdir(args.pep_folder):
print(f"PEP 文件夹不存在: {args.pep_folder}")
sys.exit(1)
for fname in os.listdir(args.pep_folder):
if fname.endswith(('.fa', '.fasta', '.faa')):
species = fname.split('.')[0]
if not species:
continue
full_path = os.path.join(args.pep_folder, fname)
if not os.path.isfile(full_path):
continue
if species in protein_path:
print(f"警告：物种名 {species} 重复，使用第一个文件。")
else:
protein_path[species] = full_path
found_species.append(species)
if not protein_path:
print(f"错误：在 {args.pep_folder} 中未找到任何蛋白文件。")
sys.exit(1)
# 3. 解析 ingroups / outgroups
focus_species_name = {}
for sp in args.ingroups.split(','):
sp = sp.strip()
if sp:
focus_species_name[sp] = True
outgroup_species_name = {}
for sp in args.outgroups.split(','):
sp = sp.strip()
if sp:
outgroup_species_name[sp] = True
# 4. 基本校验：内群/外群不能重叠，且必须存在于 PEP 文件夹的物种列表中
all_species = set(found_species)
for sp in focus_species_name:
if sp in outgroup_species_name:
print(f"错误：{sp} 同时在内群和外群中！")
sys.exit(1)
if sp not in all_species:
print(f"错误：内群物种 {sp} 不在 PEP 文件夹中。")
sys.exit(1)
for sp in outgroup_species_name:
if sp not in all_species:
print(f"错误：外群物种 {sp} 不在 PEP 文件夹中。")
sys.exit(1)
output_dir = args.output_dir
# 5. 构建精简配置字典
config = {
'tree_path': tree_path, # 只存路径，留给后续包处理
'protein_path': protein_path,
'focus_species_name': focus_species_name,
'outgroup_species_name': outgroup_species_name,
'output_dir': output_dir,
'found_species': found_species, # 可选，替代原来的 intree_check
'num_threads': args.num_threads,
}
return config
## 章老师原始的版本是设置了关注的物种，事实上，因为现在的电脑的算力资源是充分的。
## 因此完全可以考虑进行全部的blatp比对，并解析所有的分支上的结果。
def run_blastp_all_pairs(config):
"""
对 config['found_species'] 中的所有物种两两进行 BLASTP 比对。
"""
protein_path = config['protein_path']
species_list = config['found_species']
output_dir = config['output_dir']
num_threads = config['num_threads']
# 创建 BLAST 缓存目录
blast_cache_dir = os.path.join(output_dir, 'blast_cache')
os.makedirs(blast_cache_dir, exist_ok=True)
# 1. 为每个物种建库（如果 .phr 文件不存在）
for sp, fa_path in protein_path.items():
db_name = os.path.join(blast_cache_dir, sp)
# 简单的存在判断：查找 .phr 文件
if not os.path.exists(db_name + '.phr'):
print(f"[建库] {sp} <- {fa_path}")
cmd = ['makeblastdb',
'-in', fa_path,
'-dbtype', 'prot',
'-out', db_name]
try:
subprocess.run(cmd, check=True, capture_output=True, text=True)
except subprocess.CalledProcessError as e:
print(f"建库失败 {sp}:\n{e.stderr}")
raise
else:
print(f"[跳过建库] {sp} 已存在")
# 2. 两两比对
for sp_db in species_list:
db_path = os.path.join(blast_cache_dir, sp_db)
for sp_qry in species_list:
query_path = protein_path[sp_qry]
out_file = os.path.join(
blast_cache_dir,
f"{sp_qry}.2.{sp_db}_aa2aa_blast+.query"
)
print(f"[BLASTP] query={sp_qry} db={sp_db}")
cmd = [
'blastp',
'-db', db_path,
'-query', query_path,
'-out', out_file,
'-evalue', '0.00001',
'-outfmt', '7',
'-max_target_seqs', '10',
'-num_threads', str(num_threads)
]
try:
subprocess.run(cmd, check=True, capture_output=True, text=True)
except subprocess.CalledProcessError as e:
print(f"BLASTP 失败 {sp_qry} vs {sp_db}:\n{e.stderr}")
raise
def read_fasta_lengths(fasta_path):
"""读取蛋白质FASTA文件，返回 {序列ID: 长度} 字典"""
lengths = {}
for record in SeqIO.parse(fasta_path, "fasta"):
# 仅使用ID（空格前部分），对应原Perl的“space”分隔方式
seq_id = record.id.split()[0]
lengths[seq_id] = len(record.seq)
return lengths
def is_overlap(start1, stop1, start2, stop2):
"""判断两个区间 [start1,stop1] 和 [start2,stop2] 是否有重叠（均视为闭区间）"""
return max(start1, start2) <= min(stop1, stop2)
def process_blastp_results(blastp_file, query_fasta, subject_fasta,
long_copy_out, short_copy_out):
"""
解析 BLASTP 的 -outfmt 7 输出，并将命中分为长拷贝和短拷贝。
参数
----
blastp_file : str BLASTP 结果文件路径
query_fasta : str 查询物种的蛋白质FASTA文件
subject_fasta : str 目标物种的蛋白质FASTA文件
long_copy_out : str 长拷贝输出文件路径
short_copy_out : str 短拷贝输出文件路径
"""
# 获取序列长度
query_len = read_fasta_lengths(query_fasta)
subject_len = read_fasta_lengths(subject_fasta)
# 存储当前 query 的所有命中：{subject_id: (q_start, q_end, t_start, t_end)}
hits = {}
current_query = None
# 打开输出文件（句柄传给 scan_chimeric 写入）
long_fh = open(long_copy_out, 'w')
short_fh = open(short_copy_out, 'w')
try:
# 读取 BLASTP 结果
with open(blastp_file, 'r') as f:
for line in f:
# 跳过注释行
if line.startswith('#'):
# 遇到新 query 时，处理之前累积的 hits
if line.startswith('# Query: ') and current_query is not None:
scan_chimeric(current_query, hits, query_len, subject_len,
long_fh, short_fh)
hits = {}
if line.startswith('# Query: '):
current_query = line.split(': ')[1].strip()
continue
# 解析数据行
parts = line.strip().split('\t')
if len(parts) < 10:
continue
qid = parts[0]
sid = parts[1]
q_start = int(parts[6])
q_end = int(parts[7])
t_start = int(parts[8])
t_end = int(parts[9])
if qid != current_query:
# 理论上不会发生，但安全处理
if current_query is not None:
scan_chimeric(current_query, hits, query_len, subject_len,
long_fh, short_fh)
hits = {}
current_query = qid
# 合并同一 query-subject 对（取坐标的并集）
if sid in hits:
old_qs, old_qe, old_ts, old_te = hits[sid]
hits[sid] = (min(q_start, old_qs),
max(q_end, old_qe),
min(t_start, old_ts),
max(t_end, old_te))
else:
hits[sid] = (q_start, q_end, t_start, t_end)
# 处理最后一个 query
if current_query is not None and hits:
scan_chimeric(current_query, hits, query_len, subject_len,
long_fh, short_fh)
finally:
long_fh.close()
short_fh.close()
def scan_chimeric(query_id, hits, query_len_dict, subject_len_dict,long_fh, short_fh):
"""
对一个 query 的所有 subject 命中进行分析，
分出长拷贝（覆盖度>=80%）和短拷贝（去重叠后保留覆盖度>=20%的）。
"""
if query_id not in query_len_dict:
return
qlen = query_len_dict[query_id]
short_candidates = [] # 列表元素: (subject_id, q_start, q_end, t_start, t_end, coverage)
# 第一步：分类
for subject_id, (qs, qe, ts, te) in hits.items():
if subject_id not in subject_len_dict:
continue
slen = subject_len_dict[subject_id]
# 计算相对位置（归一化到0~1，保留三位小数，与原Perl一致）
qs_rel = round(qs / qlen, 3)
qe_rel = round(qe / qlen, 3)
ts_rel = round(ts / slen, 3)
te_rel = round(te / slen, 3)
coverage = round(qe_rel - qs_rel, 3)
if coverage >= 0.8:
# 长拷贝：直接输出
long_fh.write(f"{query_id}\t{subject_id}\t{qs_rel}\t{qe_rel}\t"
f"{coverage}\t{ts_rel}\t{te_rel}\n")
else:
short_candidates.append((subject_id, qs_rel, qe_rel, ts_rel, te_rel, coverage))
# 第二步：去重叠，保留覆盖度最大的短拷贝
# 按 coverage 降序排序
short_candidates.sort(key=lambda x: x[5], reverse=True)
retained = []
for item in short_candidates:
sid, qs, qe, ts, te, cov = item
# 检查是否与已保留的短拷贝重叠
overlap = False
for r in retained:
if r[0] == sid: # 同一 subject 不去重
continue
if is_overlap(qs, qe, r[1], r[2]):
overlap = True
break
if not overlap:
retained.append(item)
# 第三步：输出覆盖度 >= 0.2 的短拷贝
for sid, qs, qe, ts, te, cov in retained:
if cov >= 0.2:
short_fh.write(f"{query_id}\t{sid}\t{qs}\t{qe}\t{cov}\t{ts}\t{te}\n")
def filter_short_copies(short_copy_file, ingroup_sp, target_sp, output_dir):
"""
对短拷贝结果进行二次筛选（对应 Perl 的 "more_than_2" 逻辑）：
1. 统计每个基因的短拷贝数，过滤掉自比（query 核心名 == subject 核心名）
2. 保留有 >1 个短拷贝的基因，且确保来自不同的目标基因
3. 累计每个基因的总覆盖度，保留 >= 0.8 的
4. 去重后输出 .ids 文件
返回: (more_than_2_file, ids_file)
"""
chimeric = {} # gene_id -> count
info = {} # gene_id -> {number -> line}
with open(short_copy_file, 'r') as f:
for line in f:
line = line.strip()
if not line:
continue
parts = line.split('\t')
if len(parts) < 7:
continue
# 与 Perl 完全一致：整个行按 [\s+|.] 分割，比较索引0和索引2
svg_parts = re.split(r'[\s\+\|\.]', line)
if len(svg_parts) >= 3 and svg_parts[0] == svg_parts[2]:
continue # 跳过自比
query_id = parts[0]
if query_id not in chimeric:
chimeric[query_id] = 0
info[query_id] = {}
chimeric[query_id] += 1
info[query_id][chimeric[query_id]] = line
# 写入 more_than_2 文件（按短拷贝数降序）
out_dir_ingroup = os.path.join(output_dir, ingroup_sp)
more_than_2_file = os.path.join(
out_dir_ingroup,
f"{ingroup_sp}.2.{target_sp}.more_than_2_short_copy.ids"
)
with open(more_than_2_file, 'w') as out_f:
for gene_id, count in sorted(chimeric.items(), key=lambda x: x[1], reverse=True):
if count <= 1:
continue
# 按不同目标基因去重（与 Perl 的 tmphash 逻辑一致）
distinct = {}
for i in range(1, count + 1):
line_i = info[gene_id][i]
line_parts = line_i.split('\t')
subj_id = line_parts[1]
# 确定去重 key：如果 query_id 含 .+| 则取 subject 的第一段，否则取完整 subject_id
if '.' in line_parts[0] or '+' in line_parts[0] or '|' in line_parts[0]:
subj_key = re.split(r'[\.\+\|]', subj_id)[0]
else:
subj_key = subj_id
if subj_key not in distinct:
distinct[subj_key] = line_i
if len(distinct) > 1:
for key in sorted(distinct.keys()):
out_f.write(distinct[key] + '\n')
# 累计覆盖度，保留 >= 0.8 的基因
coverage_sum = {}
with open(more_than_2_file, 'r') as f:
for line in f:
line = line.strip()
if not line:
continue
parts = line.split('\t')
gene_id = parts[0]
cov = float(parts[4])
coverage_sum[gene_id] = coverage_sum.get(gene_id, 0.0) + cov
tmpids_file = os.path.join(out_dir_ingroup, f"{target_sp}.tmpids")
with open(tmpids_file, 'w') as f:
for gene_id, total_cov in coverage_sum.items():
if total_cov >= 0.8:
f.write(f"{gene_id}\n")
# awk '!a[$1]++' 去重
ids_file = os.path.join(out_dir_ingroup, f"{target_sp}.ids")
seen = set()
with open(tmpids_file, 'r') as fin, open(ids_file, 'w') as fout:
for line in fin:
gene_id = line.strip()
if gene_id and gene_id not in seen:
seen.add(gene_id)
fout.write(f"{gene_id}\n")
return more_than_2_file, ids_file
def process_ingroup_results(config, skip_step2=False):
"""
处理所有内群物种的 BLAST 结果，并构建 group1 / group2。
对应 Perl 中遍历 %focus_species_name 的主循环 + 第六步筛选。
参数
----
config : dict 配置字典
skip_step2 : bool 如果为 True，跳过处理，直接从已有文件构建 group1/group2
返回
-------
(group1, group2) : (list of str, list of str)
group1 格式: "target_sp:ingroup_sp/target_sp.ids"
group2 格式: "target_sp:ingroup_sp/ingroup_sp.2.target_sp.candidate_chimerics_long_copy.ids"
"""
protein_path = config['protein_path']
focus_species_name = config['focus_species_name']
output_dir = config['output_dir']
blast_cache_dir = os.path.join(output_dir, 'blast_cache')
group1 = []
group2 = []
if skip_step2:
for ingroup_sp in sorted(focus_species_name.keys()):
for target_sp in sorted(protein_path.keys()):
group1.append(f"{target_sp}:{ingroup_sp}/{target_sp}.ids")
file4 = f"{ingroup_sp}.2.{target_sp}.candidate_chimerics_long_copy.ids"
group2.append(f"{target_sp}:{ingroup_sp}/{file4}")
return group1, group2
for ingroup_sp in sorted(focus_species_name.keys()):
ingroup_dir = os.path.join(output_dir, ingroup_sp)
os.makedirs(ingroup_dir, exist_ok=True)
query_fasta = protein_path[ingroup_sp]
for target_sp in sorted(protein_path.keys()):
target_fasta = protein_path[target_sp]
blast_file = os.path.join(
blast_cache_dir,
f"{ingroup_sp}.2.{target_sp}_aa2aa_blast+.query"
)
if not os.path.isfile(blast_file):
print(f"[跳过] BLAST 结果不存在: {blast_file}")
continue
long_copy_file = os.path.join(
ingroup_dir,
f"{ingroup_sp}.2.{target_sp}.candidate_chimerics_long_copy.ids"
)
short_copy_file = os.path.join(
ingroup_dir,
f"{ingroup_sp}.2.{target_sp}.candidate_chimeric_short_copy.ids"
)
print(f"[解析 BLAST] {ingroup_sp} vs {target_sp}")
process_blastp_results(
blast_file, query_fasta, target_fasta,
long_copy_file, short_copy_file
)
# group2: 记录长拷贝文件
group2.append(
f"{target_sp}:{ingroup_sp}/"
f"{ingroup_sp}.2.{target_sp}.candidate_chimerics_long_copy.ids"
)
# 第六步：对短拷贝进行二次筛选
for target_sp in sorted(protein_path.keys()):
short_copy_file = os.path.join(
ingroup_dir,
f"{ingroup_sp}.2.{target_sp}.candidate_chimeric_short_copy.ids"
)
if not os.path.isfile(short_copy_file):
continue
print(f"[筛选短拷贝] {ingroup_sp} vs {target_sp}")
filter_short_copies(short_copy_file, ingroup_sp, target_sp, output_dir)
group1.append(f"{target_sp}:{ingroup_sp}/{target_sp}.ids")
return group1, group2
def get_tree_levels(tree, focus_species, outgroup_species):
"""
从系统发育树生成不同分类层级（对应 Perl 的 @species_in_edit_tree 迭代）。
从焦点物种向根节点逐层上移，每层产生 (ingroup, midgroup, outgroup) 三元组。
仅当 midgroup 非空时才产出（与 Perl 的 check_midgroup 逻辑一致）。
"""
all_leaves = set(tree.get_leaf_names())
outgroup_set = set(outgroup_species)
focus_leaf = tree.search_nodes(name=focus_species)[0]
node = focus_leaf
while node.up is not None:
parent = node.up
if parent.is_root():
ingroup_set = set(parent.get_leaf_names())
midgroup_set = all_leaves - ingroup_set - outgroup_set
if midgroup_set:
yield (sorted(ingroup_set), sorted(midgroup_set), sorted(outgroup_set))
# 也尝试仅焦点物种自身作为内群
ingroup_set = {focus_species}
midgroup_set = all_leaves - ingroup_set - outgroup_set
if midgroup_set:
yield (sorted(ingroup_set), sorted(midgroup_set), sorted(outgroup_set))
break
ingroup_set = set(parent.get_leaf_names())
midgroup_set = all_leaves - ingroup_set - outgroup_set
if midgroup_set:
yield (sorted(ingroup_set), sorted(midgroup_set), sorted(outgroup_set))
node = parent
def build_chimeric_matrix(config, group1, group2, focus_species):
"""
从 group1 / group2 构建嵌合基因矩阵。
返回 (all_chimeric_hash, species_chimeric_hash, candidate_orthologous)
- all_chimeric_hash: set of all candidate gene IDs
- species_chimeric_hash: {gene_id: {species: 0/1}} 短拷贝存在性
- candidate_orthologous: {gene_id: {species: 0/1}} 长拷贝存在性
"""
all_species = config['found_species']
output_dir = config['output_dir']
all_chimeric = set()
species_chimeric = {} # {gene_id: {species: 0/1}}
candidate_orthologous = {} # {gene_id: {species: 0/1}}
# 从 group1 加载短拷贝信息
for entry in group1:
species, file_name = entry.split(':', 1)
parts = file_name.split('/')
if parts[0] != focus_species:
continue
file_path = os.path.join(output_dir, file_name)
if not os.path.isfile(file_path):
continue
with open(file_path, 'r') as f:
for line in f:
gene_id = line.strip()
if not gene_id:
continue
all_chimeric.add(gene_id)
if gene_id not in species_chimeric:
species_chimeric[gene_id] = {}
species_chimeric[gene_id][species] = 1
# 初始化：所有物种的短拷贝/长拷贝都设为 0
for gene_id in all_chimeric:
for sp in all_species:
if sp not in species_chimeric[gene_id]:
species_chimeric[gene_id][sp] = 0
if gene_id not in candidate_orthologous:
candidate_orthologous[gene_id] = {}
for sp in all_species:
candidate_orthologous[gene_id][sp] = 0
# 从 group2 加载长拷贝信息
for entry in group2:
species, file_name = entry.split(':', 1)
parts = file_name.split('/')
if parts[0] != focus_species:
continue
file_path = os.path.join(output_dir, file_name)
if not os.path.isfile(file_path):
continue
with open(file_path, 'r') as f:
for line in f:
line = line.strip()
if not line:
continue
cols = line.split('\t')
gene_id = cols[0]
if gene_id in all_chimeric:
candidate_orthologous[gene_id][species] = 1
return all_chimeric, species_chimeric, candidate_orthologous
def compute_pattern(count_in_group, group_size):
"""
与 Perl 一致的 0/1 判定逻辑：
- 如果 group_size == 1: 只要有就算 1
- 如果 group_size >= 2: 需要 >= (group_size - 1) 个物种中都有才算 1
"""
if group_size <= 1:
return '1' if count_in_group >= 1 else '0'
else:
return '1' if count_in_group >= (group_size - 1) else '0'
def get_long_copy_details(gene_id, ingroup_sp, species, output_dir):
"""
从长拷贝文件中提取该基因在指定物种中的详细信息。
返回 HTML 格式的字符串（不含基因ID前缀，<br>分隔），与 Perl 的 grep+awk 一致。
"""
file_path = os.path.join(
output_dir, ingroup_sp,
f"{ingroup_sp}.2.{species}.candidate_chimerics_long_copy.ids"
)
if not os.path.isfile(file_path):
return ""
lines = []
with open(file_path, 'r') as f:
for line in f:
line = line.strip()
if line.startswith(gene_id + '\t'):
# 去掉基因ID前缀，追加 <br>（对应 Perl: $info[$i]=~s/$key//; $info[$i].="<br>"）
detail = line[len(gene_id):]
lines.append(detail + "<br>")
return ''.join(lines)
def get_short_copy_details(gene_id, ingroup_sp, species, output_dir):
"""
从 more_than_2 或原始 short_copy 文件中提取短拷贝详细信息。
返回 (details_html, lines_for_draw)。
lines_for_draw 是用于绘图的原始行列表。
"""
# 先尝试 more_than_2 文件
more_file = os.path.join(
output_dir, ingroup_sp,
f"{ingroup_sp}.2.{species}.more_than_2_short_copy.ids"
)
sc_file = os.path.join(
output_dir, ingroup_sp,
f"{ingroup_sp}.2.{species}.candidate_chimeric_short_copy.ids"
)
info_lines = []
if os.path.isfile(more_file):
with open(more_file, 'r') as f:
for line in f:
line = line.strip()
if line.startswith(gene_id + '\t'):
info_lines.append(line)
elif os.path.isfile(sc_file):
with open(sc_file, 'r') as f:
for line in f:
line = line.strip()
if line.startswith(gene_id + '\t'):
info_lines.append(line)
details = []
for line in info_lines:
detail = line[len(gene_id):]
details.append(detail + "<br>")
return ''.join(details), info_lines
def draw_chimeric_region(info_lines):
"""
对应 Perl 的 draw_chimeric_region 子程序。
根据短拷贝坐标绘制红绿相间的条形图（HTML 表格）。
"""
if not info_lines:
return ""
positions = []
for line in info_lines:
parts = line.strip().split('\t')
# 格式: query_id\tsubject_id\tqs\tqe\tcov\tts\tte
if len(parts) >= 4:
positions.append(float(parts[2])) # q_start
positions.append(float(parts[3])) # q_end
positions.sort()
html_parts = [
'<table border="0" cellspacing="0" cellpadding="0" width="100%"><tr bgcolor="red">'
]
for i, pos in enumerate(positions):
if i == 0:
percent = max(0, int(pos * 100 - 1))
else:
percent = max(0, int((pos - positions[i - 1]) * 100))
if i % 2 == 0:
# 偶数索引：空白（红色背景下的空白区域 = 非嵌合区）
html_parts.append(f'<td width="{percent}%"> </td>')
else:
# 奇数索引：绿色（嵌合区域）
html_parts.append(f'<td width="{percent}%" bgcolor="green"> </td>')
# 最后一段空白
final_percent = max(0, 100 - int(positions[-1] * 100))
html_parts.append(f'<td width="{final_percent}%"> </td>')
html_parts.append('</tr></table>')
return ''.join(html_parts)
def generate_final_output(config, group1, group2):
"""
第三步：构建嵌合基因矩阵，按树层级分类，检测最终嵌合基因，生成 HTML 报告。
对应 Perl 第 463-819 行。
"""
tree_path = config['tree_path']
focus_species_name = config['focus_species_name']
outgroup_species_name = config['outgroup_species_name']
output_dir = config['output_dir']
# 解析系统发育树
with open(tree_path, 'r') as f:
newick_str = f.read().strip()
tree = Tree(newick_str, format=1)
# 打开 all_details.txt
all_details_path = os.path.join(output_dir, 'all_details.txt')
all_details_fh = open(all_details_path, 'w')
for focus_sp in sorted(focus_species_name.keys()):
# 输出最终嵌合基因 ID 列表
final_ids_path = os.path.join(output_dir, f"{focus_sp}.final_chimeric_gene_ids")
final_ids_fh = open(final_ids_path, 'w')
# 构建矩阵
all_chimeric, species_chimeric, candidate_orthologous = \
build_chimeric_matrix(config, group1, group2, focus_sp)
# 生成 HTML
html_path = os.path.join(output_dir, f"{focus_sp}_final_chimeric_genes.html")
html_fh = open(html_path, 'w')
html_fh.write(f"""<head>
<title>Chimeric genes at {focus_sp}</title>
<meta http-equiv="content-type" content="text/html; charset=UTF-8" />
<meta name="keywords" content="chimeric gene, whole genome level scan" />
<meta name="description" content="chimeric gene scan results" />
<meta name="author" content="Chengjun Zhang" />
<meta name="copyright" content="Copyright @ Long lab" />
</head><body>
<STYLE type="text/css">
td {{white-space:nowrap}}
</STYLE>
""")
# 遍历树的不同层级
for ingroup_names, midgroup_names, outgroup_names in \
get_tree_levels(tree, focus_sp, list(outgroup_species_name.keys())):
ingroup = [x for x in ingroup_names if x]
midgroup = [x for x in midgroup_names if x]
outgroup = [x for x in outgroup_names if x]
n_ingroup = len(ingroup)
n_midgroup = len(midgroup)
n_outgroup = len(outgroup)
all_details_fh.write(f"INGROUP {' '.join(ingroup)}\n")
all_details_fh.write(f"MIDGROUP {' '.join(midgroup)}\n")
all_details_fh.write(f"OUTGROUP {' '.join(outgroup)}\n")
# HTML 表格头部
html_fh.write(f"<div><h3>{focus_sp}</h3></div>\n")
html_fh.write(
'<div><table border="1" cellspacing="1" cellpadding="1" width="100%">\n'
)
html_fh.write(
f'<TR><td>chimeric gene id</td>'
f'<td colspan={n_ingroup}>INGROUP</td>'
f'<td colspan={n_midgroup}>MIDGROUP</td>'
f'<td colspan={n_outgroup}>OUTGROUP</td></tr>\n'
f'<TR bgcolor="#F0FFFF"><td> </td>'
)
for sp in ingroup:
html_fh.write(f"<td>{sp}</td>\n")
for sp in midgroup:
html_fh.write(f"<td>{sp}</td>\n")
for sp in outgroup:
html_fh.write(f"<td>{sp}</td>\n")
html_fh.write("</tr>\n")
# 遍历每个候选嵌合基因
for gene_id in sorted(all_chimeric):
# 计算各群中的短拷贝和长拷贝计数
ingroup_sc = sum(species_chimeric[gene_id].get(sp, 0) for sp in ingroup)
midgroup_sc = sum(species_chimeric[gene_id].get(sp, 0) for sp in midgroup)
outgroup_sc = sum(species_chimeric[gene_id].get(sp, 0) for sp in outgroup)
ingroup_lc = sum(candidate_orthologous[gene_id].get(sp, 0) for sp in ingroup)
midgroup_lc = sum(candidate_orthologous[gene_id].get(sp, 0) for sp in midgroup)
outgroup_lc = sum(candidate_orthologous[gene_id].get(sp, 0) for sp in outgroup)
# 构建 0/1 模式字符串
short_pattern = (
compute_pattern(ingroup_sc, n_ingroup) +
compute_pattern(midgroup_sc, n_midgroup) +
compute_pattern(outgroup_sc, n_outgroup)
)
long_pattern = (
compute_pattern(ingroup_lc, n_ingroup) +
compute_pattern(midgroup_lc, n_midgroup) +
compute_pattern(outgroup_lc, n_outgroup)
)
# 判断是否为最终嵌合基因
is_chimeric = (
(long_pattern == "110" and short_pattern == "111") or
(long_pattern == "100" and short_pattern in ("111", "110"))
)
if not is_chimeric:
continue
final_ids_fh.write(f"{gene_id}\n")
# 输出 HTML 表格行：第一行 = 长拷贝信息
html_fh.write(f"<tr><td rowspan=2>{gene_id}</td>")
for sp in ingroup:
details = get_long_copy_details(gene_id, focus_sp, sp, output_dir)
ortho_val = candidate_orthologous[gene_id].get(sp, 0)
html_fh.write(
f'<td><a href="#">{ortho_val}</a><br>{details}</td>\n'
)
for sp in midgroup:
details = get_long_copy_details(gene_id, focus_sp, sp, output_dir)
ortho_val = candidate_orthologous[gene_id].get(sp, 0)
html_fh.write(
f'<td><a href="#">{ortho_val}</a><br>{details}</td>\n'
)
for sp in outgroup:
details = get_long_copy_details(gene_id, focus_sp, sp, output_dir)
ortho_val = candidate_orthologous[gene_id].get(sp, 0)
html_fh.write(
f'<td><a href="#">{ortho_val}</a><br>{details}</td>\n'
)
html_fh.write("</tr><tr>")
# 第二行：短拷贝信息 + 可视化
for sp in ingroup:
details, draw_lines = get_short_copy_details(
gene_id, focus_sp, sp, output_dir
)
sc_val = species_chimeric[gene_id].get(sp, 0)
if sc_val == 1:
# 短拷贝恰好 1：直接用 more_than_2 数据绘图
region_html = draw_chimeric_region(draw_lines)
html_fh.write(
f'<td><a href="#">{sc_val}</a><br>{details}\n'
f'{region_html}</td>\n'
)
else:
# 短拷贝 >1 或不存在：从原始 short_copy 文件获取绘图数据
sc_file = os.path.join(
output_dir, focus_sp,
f"{focus_sp}.2.{sp}.candidate_chimeric_short_copy.ids"
)
raw_draw_lines = []
if os.path.isfile(sc_file):
with open(sc_file, 'r') as f:
for line in f:
line = line.strip()
if line.startswith(gene_id + '\t'):
raw_draw_lines.append(line)
if raw_draw_lines:
if draw_lines:
region_html = draw_chimeric_region(draw_lines)
else:
region_html = draw_chimeric_region(raw_draw_lines[:1])
# 显示原始短拷贝信息
raw_details = []
for rl in raw_draw_lines:
raw_details.append(rl[len(gene_id):] + "<br>")
html_fh.write(
f'<td><a href="#">{sc_val}</a><br>'
f'{"".join(raw_details)}\n'
f'{region_html}</td>\n'
)
else:
html_fh.write(
f'<td><a href="#">{sc_val}</a><br>{details}</td>\n'
)
for sp in midgroup:
details, draw_lines = get_short_copy_details(
gene_id, focus_sp, sp, output_dir
)
sc_val = species_chimeric[gene_id].get(sp, 0)
if sc_val == 1:
region_html = draw_chimeric_region(draw_lines)
html_fh.write(
f'<td><a href="#">{sc_val}</a><br>{details}\n'
f'{region_html}</td>\n'
)
else:
sc_file = os.path.join(
output_dir, focus_sp,
f"{focus_sp}.2.{sp}.candidate_chimeric_short_copy.ids"
)
raw_draw_lines = []
if os.path.isfile(sc_file):
with open(sc_file, 'r') as f:
for line in f:
line = line.strip()
if line.startswith(gene_id + '\t'):
raw_draw_lines.append(line)
if raw_draw_lines:
if draw_lines:
region_html = draw_chimeric_region(draw_lines)
else:
region_html = draw_chimeric_region(raw_draw_lines[:1])
raw_details = []
for rl in raw_draw_lines:
raw_details.append(rl[len(gene_id):] + "<br>")
html_fh.write(
f'<td><a href="#">{sc_val}</a><br>'
f'{"".join(raw_details)}\n'
f'{region_html}</td>\n'
)
else:
html_fh.write(
f'<td><a href="#">{sc_val}</a><br>{details}</td>\n'
)
for sp in outgroup:
details, draw_lines = get_short_copy_details(
gene_id, focus_sp, sp, output_dir
)
sc_val = species_chimeric[gene_id].get(sp, 0)
if sc_val == 1:
region_html = draw_chimeric_region(draw_lines)
html_fh.write(
f'<td><a href="#">{sc_val}</a><br>{details}\n'
f'{region_html}</td>\n'
)
else:
sc_file = os.path.join(
output_dir, focus_sp,
f"{focus_sp}.2.{sp}.candidate_chimeric_short_copy.ids"
)
raw_draw_lines = []
if os.path.isfile(sc_file):
with open(sc_file, 'r') as f:
for line in f:
line = line.strip()
if line.startswith(gene_id + '\t'):
raw_draw_lines.append(line)
if raw_draw_lines:
if draw_lines:
region_html = draw_chimeric_region(draw_lines)
else:
region_html = draw_chimeric_region(raw_draw_lines[:1])
raw_details = []
for rl in raw_draw_lines:
raw_details.append(rl[len(gene_id):] + "<br>")
html_fh.write(
f'<td><a href="#">{sc_val}</a><br>'
f'{"".join(raw_details)}\n'
f'{region_html}</td>\n'
)
else:
html_fh.write(
f'<td><a href="#">{sc_val}</a><br>{details}</td>\n'
)
html_fh.write("</tr>\n")
html_fh.write("</table></div>\n")
html_fh.write("</body></html>\n")
html_fh.close()
final_ids_fh.close()
print(f"[完成] {focus_sp}: 最终嵌合基因 -> {final_ids_path}")
print(f"[完成] {focus_sp}: HTML 报告 -> {html_path}")
all_details_fh.close()
print(f"[完成] 全部分类明细 -> {all_details_path}")
if __name__ == "__main__":
config = parse_args_and_build_config()
print("控制台输入解析成功。")
print("树文件路径:", config['tree_path'])
print("发现物种:", config['found_species'])
print("内群:", list(config['focus_species_name'].keys()))
print("外群:", list(config['outgroup_species_name'].keys()))
print("输出目录:", config['output_dir'])
# 第一步：运行所有两两 BLASTP 比对
run_blastp_all_pairs(config)
# 第二步：解析 BLAST 结果并筛选短拷贝
group1, group2 = process_ingroup_results(config)
print(f"\n处理完成。group1: {len(group1)} 条, group2: {len(group2)} 条")
# 第三步：树层级分类 + 嵌合基因检测 + HTML 报告
generate_final_output(config, group1, group2)我跑的out忘记指定了，明天重新跑吧，现在来看blastp消耗了大量的时间，还是很有必要处理一下的。
对blatsp版本不通的解决办法
在数据集相同的情况下，让python版本的程序和perl版本的程序调用相同的blastp文件。
blastp命令报错
我也遇到了这个错误，上github查看，果然有:
https://github.com/zhangcj2022/GriffinDetector/issues/4
这个是因为2.2.31版本以及我们目前最新版本的blastp建库命令已经不是这个命令
$cmd="$makeblastdb -in $dbfile -hash_index -dbtype prot -parse_seqids -gi_mask";手动修改成这个命令
makeblastdb -in $dbfile -dbtype prot还有一个bug，他这里的配置文件的文件名必须是锁死的GriffinDetector.ctl
但是命令行调用的时候又要手动输入，所以很容易让人以为这个文件是可以随意修改的。
comm -23 <(sort -u /project1/data/identified_chimerics/grandisfinal_chimeric_gene_ids) <(sort -u /project1/data/used_for_fusion/output/grandis.final_chimeric_gene_ids)Perl版本与python版本结果的匹配问题
perl版本比最后一个python版本多了4个基因(去冗余后)，先检查python侧
for G in evm_model_PTG001531L_56 evm_model_PTG001984L_93 evm_model_PTG004650L_4 evm_model_PTG011191L_11; do
echo "=== $G ==="
for f in grandis.2.*.candidate_chimeric_short_copy.ids; do
count=$(grep -c "$G" "$f" 2>/dev/null)
[ "$count" -gt 0 ] && echo " $f: $count 行"
done
echo "--- long_copy ---"
grep -l "$G" grandis.2.*.candidate_chimerics_long_copy.ids 2>/dev/null || echo " 无"
echo "--- more_than_2 ---"
grep -l "$G" grandis.2.*.more_than_2_short_copy.ids 2>/dev/null || echo " 无"
echo "--- .ids ---"
grep -l "$G" *.ids 2>/dev/null | grep -v candidate | grep -v more | grep -v tmp || echo " 无"
echo
done
检测perl一侧的输出
(gffrin) root@53d03939188f:/project1/data/identified_chimerics/grandis# for G in evm_model_PTG001531L_56 evm_model_PTG004650L_4 evm_model_PTG011191L_11; do
echo "=== $G ==="
for f in grandis2*candidate_chimeric_short_copy.ids; do
count=$(grep -c "$G" "$f" 2>/dev/null)
[ "$count" -gt 0 ] && echo " $f: $count 行"
done
echo
done
=== evm_model_PTG001531L_56 ===
grandis2californicacandidate_chimeric_short_copy.ids: 2 行
grandis2chinese_pinecandidate_chimeric_short_copy.ids: 2 行
grandis2fargesiicandidate_chimeric_short_copy.ids: 2 行
grandis2grandiscandidate_chimeric_short_copy.ids: 2 行
grandis2jackiicandidate_chimeric_short_copy.ids: 2 行
grandis2nuciferacandidate_chimeric_short_copy.ids: 2 行
grandis2taxifoliacandidate_chimeric_short_copy.ids: 2 行
grandis2taxus_chinensiscandidate_chimeric_short_copy.ids: 2 行
=== evm_model_PTG004650L_4 ===
grandis2californicacandidate_chimeric_short_copy.ids: 2 行
grandis2chinese_pinecandidate_chimeric_short_copy.ids: 2 行
grandis2fargesiicandidate_chimeric_short_copy.ids: 2 行
grandis2grandiscandidate_chimeric_short_copy.ids: 2 行
grandis2jackiicandidate_chimeric_short_copy.ids: 2 行
grandis2nuciferacandidate_chimeric_short_copy.ids: 2 行
grandis2taxifoliacandidate_chimeric_short_copy.ids: 2 行
grandis2taxus_chinensiscandidate_chimeric_short_copy.ids: 2 行
=== evm_model_PTG011191L_11 ===
grandis2californicacandidate_chimeric_short_copy.ids: 2 行
grandis2chinese_pinecandidate_chimeric_short_copy.ids: 2 行
grandis2fargesiicandidate_chimeric_short_copy.ids: 2 行
grandis2grandiscandidate_chimeric_short_copy.ids: 2 行
grandis2jackiicandidate_chimeric_short_copy.ids: 2 行
grandis2nuciferacandidate_chimeric_short_copy.ids: 2 行
grandis2taxifoliacandidate_chimeric_short_copy.ids: 2 行
grandis2taxus_chinensiscandidate_chimeric_short_copy.ids: 2 行(gffrin) root@53d03939188f:/project1/data/identified_chimerics/grandis# for G in evm_model_PTG001531L_56 evm_model_PTG001984L_93 evm_model_PTG004650L_4 evm_model_PTG011191L_11; do
echo "=== $G (californica) ==="
grep "$G" grandis2californicacandidate_chimeric_short_copy.ids
echo
done
=== evm_model_PTG001531L_56 (californica) ===
evm_model_PTG001531L_56 californica_JS2-8gene32672 0.009 0.486 0.477 0.961 0.997
evm_model_PTG001531L_56 californica_JS2-8gene886 0.009 0.486 0.477 0.788 0.981
=== evm_model_PTG001984L_93 (californica) ===
evm_model_PTG001984L_93 californica_JS2-6gene27489 0.007 0.472 0.465 0.505 0.996
evm_model_PTG001984L_93 californica_JS2-8gene17272 0.57 0.972 0.402 0.119 0.661
evm_model_PTG001984L_93 californica_JS2-6gene27489 0.007 0.472 0.465 0.505 0.996
evm_model_PTG001984L_93 californica_JS2-8gene17272 0.57 0.972 0.402 0.119 0.661
=== evm_model_PTG004650L_4 (californica) ===
evm_model_PTG004650L_4 californica_JS2-6gene22450 0.138 0.713 0.575 0.554 0.76
evm_model_PTG004650L_4 californica_JS2-6gene22450 0.138 0.713 0.575 0.554 0.76
=== evm_model_PTG011191L_11 (californica) ===
evm_model_PTG011191L_11 californica_JS2-8gene9117 0.02 0.56 0.54 0.955 0.971
evm_model_PTG011191L_11 californica_JS2-8gene9115 0.02 0.56 0.54 0.768 0.851本质上就是一个query比对到两个id不同的，但是比对区域完全相同的基因，对不对？这种情况理论上应该是要删除的，但是在perl版本里被保留了。你能找到perl版本对应的代码，并告诉我为什么吗？
一直到最后都没能复刻这四个基因，但是可以很明确的是，这4个基因是假阳性，没有复刻的必要。
Iqtree 所构建进化树
这个树的和下面这个树的区别是买麻藤和百岁兰分在中间，银杏在外面，两个被子植物中有一个和裸子植物聚类到一起了。
群体数据集的可用性：
为了将来计算群体的选择压力，我想要搞清楚，这些数据集里面，有群体vcf的数据集有哪些。
很多物种在文章中做了群体工作，但是未必有vcf数据集。
香榧肯定是有的，Tgrandis；
Ortherfinder所构建进化树
这两个树都是可用的，
更新计划
- 考虑短拷贝之间的同源性，用于判断fusion事件是否是独立起源的。
- 分branchers遍历所有的fusion事件，构建矩阵，这样一来，如果不同物种的fusion基因内部考虑同源性，那么就可以把基因fusion事件在不同branch上的发生事件和特定物种上的丢失事件找出来。
- 根据gff找到fusion在不同染色体上的发生位置，并判定其发生方式(染色体重排？逆转录转座？)。
原始的fusion代码考虑了基因在不同染色体上的位置，但是仅仅用作提取注释之用；
22年gb的文章仅仅做了理论分类，而没有具体分析起源的发生机制。
由于裸子植物本身保持着良好的基因组共线性，在对染色体臂进行合理编号的情况下，或许可以对发生的机制进行分析。
短拷贝保守性检验
一个真实的融合基因 = 多个祖先基因片段拼接而成。每个片段在外群物种中应有保守的直系同源。短拷贝在这些外群物种中的位置应当是聚集的，而非散乱分布。
同一个 branch 上的融合基因
│
▼
按短拷贝位置聚类 → 得到每个基因的"位置签名"
│
▼
相同签名的基因 ──→ 同一祖先融合事件（各物种继承或独立dup）
不同签名的基因 ──→ 独立起源的融合事件
进化树的动态划分
原始的划分方法是这样的：
对，是动态划分的。没有阈值——完全靠树的拓扑结构逐层压缩。过程如下：
原始树（简化版，假设 focus=Tchinensis，外群=Mpolymorpha）：
第 1 轮：regex 找到(Tchinensis,Sgiganteum)→ 标为_a1
树变成：(Mpolymorpha,(Atrichopoda,(Gmontanum,(Gbiloba,(Ptabuliformis,_a1)))))
划分结果：
ingroup：Tchinensis（focus 本身）
midgroup：Sgiganteum、Ptabuliformis、Gbiloba、Gmontanum、Atrichopoda（树上夹在 focus 和外群之间的所有物种）
outgroup：Mpolymorpha
第 2 轮：regex 找到(Ptabuliformis,_a1)→ 标为_b1
划分：
ingroup：Tchinensis、Sgiganteum、Ptabuliformis（_b1展开后包含 focus）
midgroup：Gbiloba、Gmontanum、Atrichopoda
outgroup：Mpolymorpha
第 3 轮：(Gbiloba,_b1)→_c1，ingroup 继续扩大...
直到 midgroup 为空（所有物种都进了 ingroup）。
核心规律：
ingroup 每次都向根方向"吞掉"一个物种
midgroup 逐渐缩小
同一基因在任意一层匹配 pattern 就算嵌合
这就是为什么 midgroup 越大 → 融合越古老——那时候 ingroup 还很小，外群之外的物种都还在中群里。
但是啊，我觉得这样的划分方法有问题
但是啊，这里有个问题，有时候，对于一个内fucusion物种，midgroup中某两个物种可能是等同的，即处于系统发育树的相同拓扑位置，所以我开始觉得原始perl代码这样划分是有失平衡的。
我觉得应当按照相对于fucus物种的节点来划分，你觉得对吗?
对了，顺便看一下原始perl代码有什么问题：
现在版本的代码的处理方式，以下面这个树为例：
[v8 tree] focus=Tgrandis 祖先层级数=7
Level 1: +4 species ingroup=(4) | midgroup=(15)['Cgigantea', 'Cjaponica', 'Clanceolata']...
DEBUG: [Tgrandis] level ingroup=['Pseudotaxus_chienii', 'Tchinensis', 'Tgrandis', 'Twallichiana'] midgroup=['Cgigantea', 'Cjaponica', 'Clanceolata', 'Cpanzhihuaensis', 'Fhodginsii', 'Gbiloba', 'Gmontanum', 'Mglyptostroboides', 'Pdensiflora', 'Platycladus_orientalis', 'Pmassoniana', 'Ptabuliformis', 'Sgiganteum', 'Tdistichum', 'Wmirabilis'] -> total_chimeric=5403
DEBUG: [Tgrandis] gene=TgF01G00100-mRNA1 sc=(0/4,1/15,0/3) lc=(1/4,1/15,0/3) -> short=000 long=000 FAIL
DEBUG: [Tgrandis] gene=TgF01G00240-mRNA1 sc=(0/4,1/15,0/3) lc=(3/4,14/15,3/3) -> short=000 long=011 FAIL
DEBUG: [Tgrandis] gene=TgF01G00260-mRNA1 sc=(0/4,1/15,0/3) lc=(3/4,15/15,3/3) -> short=000 long=011 FAIL
DEBUG: [Tgrandis] level ingroup=['Pseudotaxus_chienii', 'Tchinensis', 'Tgrandis', 'Twallichiana'] -> pass/total = 0/5403
Level 2: +8 species ingroup=(12) | midgroup=(7)['Cpanzhihuaensis', 'Gbiloba', 'Gmontanum']...
DEBUG: [Tgrandis] level ingroup=['Cgigantea', 'Cjaponica', 'Clanceolata', 'Fhodginsii', 'Mglyptostroboides', 'Platycladus_orientalis', 'Pseudotaxus_chienii', 'Sgiganteum', 'Tchinensis', 'Tdistichum', 'Tgrandis', 'Twallichiana'] midgroup=['Cpanzhihuaensis', 'Gbiloba', 'Gmontanum', 'Pdensiflora', 'Pmassoniana', 'Ptabuliformis', 'Wmirabilis'] -> total_chimeric=5403
DEBUG: [Tgrandis] gene=TgF01G00100-mRNA1 sc=(1/12,0/7,0/3) lc=(1/12,1/7,0/3) -> short=000 long=000 FAIL
DEBUG: [Tgrandis] gene=TgF01G00240-mRNA1 sc=(0/12,1/7,0/3) lc=(10/12,7/7,3/3) -> short=000 long=111 FAIL
DEBUG: [Tgrandis] gene=TgF01G00260-mRNA1 sc=(1/12,0/7,0/3) lc=(11/12,7/7,3/3) -> short=000 long=111 FAIL
DEBUG: [Tgrandis] level ingroup=['Cgigantea', 'Cjaponica', 'Clanceolata', 'Fhodginsii', 'Mglyptostroboides', 'Platycladus_orientalis', 'Pseudotaxus_chienii', 'Sgiganteum', 'Tchinensis', 'Tdistichum', 'Tgrandis', 'Twallichiana'] -> pass/total = 0/5403
Level 3: +3 species ingroup=(15) | midgroup=(4)['Cpanzhihuaensis', 'Gbiloba', 'Gmontanum']...
DEBUG: [Tgrandis] level ingroup=['Cgigantea', 'Cjaponica', 'Clanceolata', 'Fhodginsii', 'Mglyptostroboides', 'Pdensiflora', 'Platycladus_orientalis', 'Pmassoniana', 'Pseudotaxus_chienii', 'Ptabuliformis', 'Sgiganteum', 'Tchinensis', 'Tdistichum', 'Tgrandis', 'Twallichiana'] midgroup=['Cpanzhihuaensis', 'Gbiloba', 'Gmontanum', 'Wmirabilis'] -> total_chimeric=5403
DEBUG: [Tgrandis] gene=TgF01G00100-mRNA1 sc=(1/15,0/4,0/3) lc=(1/15,1/4,0/3) -> short=000 long=000 FAIL
DEBUG: [Tgrandis] gene=TgF01G00240-mRNA1 sc=(1/15,0/4,0/3) lc=(13/15,4/4,3/3) -> short=000 long=111 FAIL
DEBUG: [Tgrandis] gene=TgF01G00260-mRNA1 sc=(1/15,0/4,0/3) lc=(14/15,4/4,3/3) -> short=000 long=111 FAIL
DEBUG: [Tgrandis] level ingroup=['Cgigantea', 'Cjaponica', 'Clanceolata', 'Fhodginsii', 'Mglyptostroboides', 'Pdensiflora', 'Platycladus_orientalis', 'Pmassoniana', 'Pseudotaxus_chienii', 'Ptabuliformis', 'Sgiganteum', 'Tchinensis', 'Tdistichum', 'Tgrandis', 'Twallichiana'] -> pass/total = 0/5403
Level 4: +2 species ingroup=(17) | midgroup=(2)['Gmontanum', 'Wmirabilis']
DEBUG: [Tgrandis] level ingroup=['Cgigantea', 'Cjaponica', 'Clanceolata', 'Cpanzhihuaensis', 'Fhodginsii', 'Gbiloba', 'Mglyptostroboides', 'Pdensiflora', 'Platycladus_orientalis', 'Pmassoniana', 'Pseudotaxus_chienii', 'Ptabuliformis', 'Sgiganteum', 'Tchinensis', 'Tdistichum', 'Tgrandis', 'Twallichiana'] midgroup=['Gmontanum', 'Wmirabilis'] -> total_chimeric=5403
DEBUG: [Tgrandis] gene=TgF01G00100-mRNA1 sc=(1/17,0/2,0/3) lc=(2/17,0/2,0/3) -> short=000 long=000 FAIL
DEBUG: [Tgrandis] gene=TgF01G00240-mRNA1 sc=(1/17,0/2,0/3) lc=(15/17,2/2,3/3) -> short=000 long=111 FAIL
DEBUG: [Tgrandis] gene=TgF01G00260-mRNA1 sc=(1/17,0/2,0/3) lc=(16/17,2/2,3/3) -> short=000 long=111 FAIL
DEBUG: [Tgrandis] level ingroup=['Cgigantea', 'Cjaponica', 'Clanceolata', 'Cpanzhihuaensis', 'Fhodginsii', 'Gbiloba', 'Mglyptostroboides', 'Pdensiflora', 'Platycladus_orientalis', 'Pmassoniana', 'Pseudotaxus_chienii', 'Ptabuliformis', 'Sgiganteum', 'Tchinensis', 'Tdistichum', 'Tgrandis', 'Twallichiana'] -> pass/total = 0/5403在v8版本的基础上我又想出了一个新的版本：
即双扩散迭代和S层。特别是S层，被专门用于寻找那些物种特异的fusion事件。
这个迭代虽然能检测出物种特异性分质上的fusion事件。
但是，对于fusion发生branch和发生类型（短拷贝保留or不保留）的处理上，带来了困扰。
但是，现在我又有了更好的办法：
是的，我们能不能用两层划分法？
只划分内群和外群？
然后在分析的时候，可以放弃一些物种。
具体来说，是这样？
以Tgrandis为例，第一步我们只考虑最近外类群，即第一个祖先节点上除去Tgrandis后剩下的三个物种。
然后把这三个物种作为外群。
如果long=10，short=01，则代表fusion发生在最近谱系上，也就是说物种特异的，同时短拷贝丢失或者消耗。
如果long=10，short=11，则代表fusion发生在在最近谱系上，同时短拷贝保留。
如果long=11 short=11，则代表fusion可能发生在更古老的分支上。
这个时候需要向外延申，再把Tgrandis+最近外类群，也就是这四个物种作为内群。
把外侧第一个近源group作为外群。
你理解我的思路了吗？
这样一来，我们能更轻松的对与fusion事件进行brchan划分和短拷贝存续统计，对不对。
这个框架的首要优点是，他是一个很美观的框架，任何一个清晰的东西，都应当是简洁而美观的。
- 不考虑物种特异性丢失的情况下，任何一个fusion事件都应当发生在特定的branch上，而这个branch内部，应当保留long_copy，branch外部，没有longe_copy
- 而这个brach以外的外类群，应当保留short_copy。
- 如果这个branch内部没有short copy呢？那就是说明short_copy丢失了，或者在融合阶段就被消耗掉了。
所以：两群融合判定的核心原则
1. 分支定位原则
任何融合事件都发生在系统发育树的特定分支上。该分支内部的物种保留长拷贝（long copy），分支外部的物种不保留长拷贝。因此 long = "10"（内群有、外群无）是融合发生在当前分支的充要条件。
2. 外群证据原则
融合事件的亲本基因早于融合分支存在。分支外部的近缘物种（外群）应当保留该亲本基因的短拷贝碎片（short copy），作为追溯融合事件的进化证据。
3. 短拷贝丢失原则
如果融合分支内部（内群）没有检出短拷贝，说明亲本基因的碎片在融合发生时已被消耗（如通过重组、缺失等机制），或在融合后的进化过程中已完全丢失。这对应 short = "01"——短拷贝仅在外群保留，内群已无。
这三条直接翻译成代码逻辑就是：
long="10" short="11" → 融合在此分支, 亲本碎片仍存 (retained)
long="10" short="01" → 融合在此分支, 亲本碎片已丢 (lost)
long="11" short="11" → 融合早于此分支, 扩大内群继续追溯所以只有到这一步为止，这个代码的逻辑判定才会是准确的。
- 这个代码完美融合了所有外群，内群，和中群。
- 以上框架，还有利于我们对一个分支内部的，物种间的，谱系特异性的long_copy进行合并。
- 还有利于我们对谱系内部的long_copy再次复制的事件，进行更深入的分析。
此外，我还发现：
- 我发现对于long_copy=1的标准要放的很宽松，只要有一个物种，就视为1，只有全部物种都没有，才能计算为0.
- 原始perl代码是怎么处理的。
所以long_copy应当使用宽松阈值，short_copy应当使用严格阈值。
Fusion发生branch标记与整合
这个部分具体是指：在获得fusion基因的融合时间后，对不同物种的，同一fusion事件产生的直系同源基因进行整合，从而标注出不同branch上发生的fusion事件的总数。
基因id处理逻辑
这个很简单，就是biopython取代原始章老师的错误划分方式，他那个划分方式对于部分物种来说，会有bug。
对fusion后fusion产生基因复制事件进行标记（还没做）
要不要加入ks标定呢？加入ks标定有一个好处，可以对每个融合基因的复制时间进行标定。
阶段性总结
以上内容包含了Fusion软件的全部，榧属融合基因的鉴定部分，接下来的内容主要是为了达到一个文章发表框架而填充内容所做的，具体而言。整个工作可以划分为以下几个部分。
- Fusion工具的升级
- 裸子植物的全部fusion事件鉴定
- fusion基因的蛋白质结构（含IDR）
蛋白质结构不急于跑。我们先看4
- fusion基因的群体动态选择压力
1. 融合基因在多个物种间（L1+ 共享事件）
可行性：低。只有少数 L1+ 共享融合事件有跨物种直系同源——每个 clade 几个到几十个。
能算的：每对物种的 pairwise Ka/Ks
dN/dS 含义
< 1 融合蛋白在物种分化后一直被保守，功能稳定
≈ 1 纯中性漂变，功能不受约束（在退化？）
> 1 某条 lineage 上融合蛋白加速演化
局限：dN/dS 是两个物种的平均值，无法区分"是在 A 变快了还是 B 变慢了"。
2. 融合基因 vs 短拷贝（同物种内，旁系同源）
可行性：高。每个融合基因都有 short copy（亲本），全在同物种内，不需跨物种。
能算的：每个融合基因 vs 亲本 A、vs 亲本 B
dN/dS 含义
< 1 融合后蛋白未退化，受净化选择约束 — 基因仍然有功能
≈ 1 融合后蛋白在漂变 — 可能是假基因化过程中
> 1 融合后蛋白被正选择推向新功能 — 极其重要的信号
这个能做批量：28 个融合基因 × 2 个亲本 = 56 对 pairwise dN/dS，一网打尽。
3. 群体水平正选择
可行性：需要多个体 SNP 数据（VCF 或重测序 reads）
能算的：
指标 含义
Tajima's D < 0 群体中稀有等位基因过剩 → 近期净化选择或选择性清除
融合基因区域的 π 显著低于基因组背景 融合基因在群体中被保守
McDonald-Kreitman 检验 融合基因在种内多态 vs 种间固定的比例偏离中性预期
没有群体数据的话，这个做不了。对于短拷贝重叠逻辑的修改
原始版本的代码中，对于短拷贝重叠情况，章老师使用了有重叠就删除较短者的思路。
在这版本的代码中，我发现了这个问题，我对于短拷贝做出了调整，我要求重叠达到70%才能删除一个短拷贝，在修改之后，我们找到了三个