基于Python统计多FASTA文件中不同分组的Busco基因数量

先梳理清楚需求：我们手里有一个多FASTA格式的序列文件，其中Busco基因的ID以>EOG开头，需要完成去重统计（同一Busco基因不管后缀编号如何只算一次），还要按两个维度拆分统计：

• 按物种编号：统计ID第一个数字段为_0035、_0042的Busco基因数量
• 按ID数字组合：分别统计_0035_0042、_0042_0035、_0042_0042、_0035_0035这四组的Busco基因数量

原脚本的优化方向 #

你原来的脚本已经实现了基础的Busco基因去重计数，但缺少维度拆分的逻辑。我们可以基于原思路扩展，用字典跟踪每个Busco基因对应的所有数字组合，这样就能轻松实现多维度统计。

改进后的Python脚本 #

from collections import defaultdict

# 存储每个Busco基因对应的所有数字组合（用集合避免重复）
busco_groups = defaultdict(set)

with open("concatenate_with_busco_names_0035_0042_aa.fa") as f:
    for line in f:
        line = line.strip()
        if line.startswith('>EOG'):
            # 拆分ID，提取核心信息
            parts = line.split('_')
            busco_name = parts[0][4:]  # 去掉">EOG"前缀
            num_pair = f"_{parts[1]}_{parts[2]}"  # 提取数字组合，比如_0035_0042
            busco_groups[busco_name].add(num_pair)

# 初始化统计变量
total_busco = len(busco_groups)
specie_0035 = 0
specie_0042 = 0
pair_0035_0042 = 0
pair_0042_0035 = 0
pair_0042_0042 = 0
pair_0035_0035 = 0

# 遍历每个Busco基因的组合集合，完成多维度统计
for busco, pairs in busco_groups.items():
    # 统计物种维度：只要有一个组合以对应编号开头，就计入该物种
    if any(p.startswith('_0035') for p in pairs):
        specie_0035 += 1
    if any(p.startswith('_0042') for p in pairs):
        specie_0042 += 1
    # 统计数字组合维度：只要该Busco出现在对应组合中，就计入统计
    if '_0035_0042' in pairs:
        pair_0035_0042 += 1
    if '_0042_0035' in pairs:
        pair_0042_0035 += 1
    if '_0042_0042' in pairs:
        pair_0042_0042 += 1
    if '_0035_0035' in pairs:
        pair_0035_0035 += 1

# 输出最终统计结果
print(f"Total busco: {total_busco} (I count only once if the >busco is present even if _number are different)")
print(f"Total busco for the specie _0035 (_0035_0042 and _0035_0035) : {specie_0035}")
print(f"Total busco for the specie _0042 (_0042_0042 and _0042_0035) : {specie_0042}")
print(f"Total busco for the specific specie _0035_0042 : {pair_0035_0042}")
print(f"Total busco for the specific specie _0042_0035 : {pair_0042_0035}")
print(f"Total busco for the specific specie _0042_0042 : {pair_0042_0042}")
print(f"Total busco for the specific specie _0035_0035 : {pair_0035_0035}")

脚本逻辑说明 #

1. 用defaultdict(set)存储每个Busco基因的数字组合，确保同一组合不会重复记录
2. 遍历FASTA文件时，只处理Busco基因的ID行，拆分提取关键信息
3. 统计阶段：
   • 物种维度：只要该Busco基因有一个组合以目标编号开头，就计入对应物种的统计
   • 数字组合维度：只要该Busco基因出现在对应组合中，就计入该组合的统计

示例输出（匹配你给出的预期结果） #

Total busco: 5 (I count only once if the >busco is present even if _number are different)
Total busco for the specie _0035 (_0035_0042 and _0035_0035) : 3
Total busco for the specie _0042 (_0042_0042 and _0042_0035) : 4
Total busco for the specific specie _0035_0042 : 3
Total busco for the specific specie _0042_0035 : 0
Total busco for the specific specie _0042_0042 : 4
Total busco for the specific specie _0035_0035 : 0

内容的提问来源于stack exchange，提问作者Grendel