筛选同一终止密码子下最长ORF的Python脚本问题求助

我有一组直系同源序列文件，其中多个蛋白带有相同的终止密码子，需要筛选出同一菌株、同一终止密码子下的最长ORF。编程经验有限，尝试用嵌套字典（层级：all_orfs = {菌株1: {终止密码子1: {长度1: 标识符1, 长度2: 标识符2}}}）实现，但未成功。比如需筛选菌株NZ_AP022573.1_Mycobacterium_saskatchewanense终止密码子3965415下的最长ORF（ORF_150482，长度89）。

原脚本代码 #

def o_smorf_gatherer(folder, outdir):
    files = os.listdir(folder)

    for file in files:
        all_orfs = {}
        records = SeqIO.parse(f'{folder}/{file}', 'fasta')
        for record in records:
            print(file)
            # 原脚本错误：未定义identifier，应该用record.id
            identifier_parts = identifier.split('_')
            # print(identifier_parts)
            species_strain = '_'.join(identifier_parts[-4:])
            # print(species_strain)
            coordinates = identifier_parts[5]
            codon = coordinates.split('-')
            # print(codon)
            stop_codon = codon[-1]
            # print(stop_codon)
            length = identifier_parts[3]
            # print(length)
            if species_strain not in all_orfs:
                all_orfs[species_strain] = {}
            if stop_codon not in all_orfs[species_strain]:
                all_orfs[species_strain][stop_codon] = {}
            if length not in all_orfs[species_strain][stop_codon]:
                all_orfs[species_strain][stop_codon][length] = {}
            if identifier not in all_orfs[species_strain][stop_codon][length]:
                all_orfs[species_strain][stop_codon][length] = [identifier]
            print(all_orfs)

            for stop_codon in all_orfs[species_strain]:
                longest_orf = None
                longest_seq = 0
                sequence_length = all_orfs[species_strain][stop_codon]
            #     if sequence_length > longest_seq:
            #         longest_seq = sequence_length
            #         longest_orf = all_orfs[species_strain][stop_codon][length]
            #         # print(longest_orf)


            # Now I need to access the length in the dictio
            # print(all_orfs[species_strain][stop_codon])
            biggest_orf = max(all_orfs[species_strain][stop_codon])
            print(biggest_orf)
            unique_orfs = {}
            if biggest_orf not in unique_orfs:
                unique_orfs[biggest_orf] = []
                print(unique_orfs)

FASTA输入示例 #

> ORF_150486_LEN_77_COORD_3965185-3965415_forward_NZ_AP022573.1_Mycobacterium_saskatchewanense
> ORF_150484_LEN_81_COORD_3965173-3965415_forward_NZ_AP022573.1_Mycobacterium_saskatchewanense
> ORF_150483_LEN_84_COORD_3965164-3965415_forward_NZ_AP022573.1_Mycobacterium_saskatchewanense
> ORF_150482_LEN_89_COORD_3965149-3965415_forward_NZ_AP022573.1_Mycobacterium_saskatchewanense
> ORF_260114_LEN_92_COORD_5871599-5871874_forward_NZ_CP070348.1_Mycolicibacterium_boenickei
> ORF_260116_LEN_60_COORD_5871695-5871874_forward_NZ_CP070348.1_Mycolicibacterium_boenickei
> ORF_120806_LEN_59_COORD_2912055-2911879_reverse_NZ_AP022612.1_Mycolicibacterium_confluentis
> ORF_120805_LEN_61_COORD_2912061-2911879_reverse_NZ_AP022612.1_Mycolicibacterium_confluentis
> ORF_120803_LEN_69_COORD_2912085-2911879_reverse_NZ_AP022612.1_Mycolicibacterium_confluentis
> ORF_118543_LEN_98_COORD_3460732-3460439_reverse_NZ_AP022560.1_Mycolicibacterium_moriokaense
> ORF_112065_LEN_55_COORD_2618565-2618729_forward_NZ_AP022607.1_Mycobacterium_branderi
> ORF_112067_LEN_53_COORD_2618571-2618729_forward_NZ_AP022607.1_Mycobacterium_branderi
> ORF_112064_LEN_57_COORD_2618559-2618729_forward_NZ_AP022607.1_Mycobacterium_branderi
> ORF_112062_LEN_65_COORD_2618535-2618729_forward_NZ_AP022607.1_Mycobacterium_branderi
> ORF_224113_LEN_59_COORD_900663-900487_reverse_NZ_AP022608.1_Mycolicibacterium_gadium
> ORF_224110_LEN_67_COORD_900687-900487_reverse_NZ_AP022608.1_Mycolicibacterium_gadium
> ORF_79139_LEN_59_COORD_1475341-1475517_forward_NZ_CP007220.1_Mycobacteroides_chelonae
> ORF_116852_LEN_77_COORD_3581529-3581299_reverse_NZ_AP022587.1_Mycobacterium_stomatepiae
> ORF_116851_LEN_79_COORD_3581535-3581299_reverse_NZ_AP022587.1_Mycobacterium_stomatepiae
> ORF_116850_LEN_81_COORD_3581541-3581299_reverse_NZ_AP022587.1_Mycobacterium_stomatepiae
> ORF_116849_LEN_84_COORD_3581550-3581299_reverse_NZ_AP022587.1_Mycobacterium_stomatepiae
> ORF_178001_LEN_86_COORD_4676508-4676251_reverse_NZ_CP078145.1_Nocardia_iowensis
> ORF_178000_LEN_89_COORD_4676517-4676251_reverse_NZ_CP078145.1_Nocardia_iowensis
> ORF_219647_LEN_84_COORD_621170-620919_reverse_NZ_AP022616.1_Mycolicibacterium_phocaicum
> ORF_219646_LEN_86_COORD_621176-620919_reverse_NZ_AP022616.1_Mycolicibacterium_phocaicum
> ORF_219645_LEN_89_COORD_621185-620919_reverse_NZ_AP022616.1_Mycolicibacterium_phocaicum
> ORF_219644_LEN_94_COORD_621200-620919_reverse_NZ_AP022616.1_Mycolicibacterium_phocaicum
> ORF_219649_LEN_54_COORD_621080-620919_reverse_NZ_AP022616.1_Mycolicibacterium_phocaicum
> ORF_40651_LEN_77_COORD_5220442-5220212_reverse_NZ_CP089224.1_Mycobacterium_ostraviense
> ORF_40650_LEN_79_COORD_5220448-5220212_reverse_NZ_CP089224.1_Mycobacterium_ostraviense
> ORF_40646_LEN_89_COORD_5220478-5220212_reverse_NZ_CP089224.1_Mycobacterium_ostraviense
> ORF_40648_LEN_84_COORD_5220463-5220212_reverse_NZ_CP089224.1_Mycobacterium_ostraviense
> ORF_40649_LEN_81_COORD_5220454-5220212_reverse_NZ_CP089224.1_Mycobacterium_ostraviense
> ORF_26333_LEN_47_COORD_546668-546808_forward_NZ_CP059165.1_Mycobacterium_vicinigordonae

修正后的脚本及说明 #

原脚本存在几个关键问题：未正确获取FASTA记录的标识符、长度以字符串存储无法正确比较、逻辑顺序错误（处理单个记录就计算最长）。以下是修正后的版本：

import os
from Bio import SeqIO

def o_smorf_gatherer(folder, outdir):
    # 创建输出目录（如果不存在）
    if not os.path.exists(outdir):
        os.makedirs(outdir)
    
    files = os.listdir(folder)
    for file in files:
        # 初始化存储结构：{菌株: {终止密码子: [(长度, 记录), ...]}}
        all_orfs = {}
        file_path = os.path.join(folder, file)
        records = SeqIO.parse(file_path, 'fasta')
        
        for record in records:
            # 处理FASTA标题，去除开头空格后分割
            identifier = record.id.strip()
            identifier_parts = identifier.split('_')
            
            # 提取菌株信息（最后4个部分）
            species_strain = '_'.join(identifier_parts[-4:])
            # 提取终止密码子坐标
            coord_part = identifier_parts[5]
            stop_codon = coord_part.split('-')[-1]
            # 提取长度并转为整数
            length = int(identifier_parts[3])
            
            # 填充字典
            if species_strain not in all_orfs:
                all_orfs[species_strain] = {}
            if stop_codon not in all_orfs[species_strain]:
                all_orfs[species_strain][stop_codon] = []
            # 将长度和记录存入列表
            all_orfs[species_strain][stop_codon].append((length, record))
        
        # 筛选每个菌株-终止密码子组的最长ORF
        output_file = os.path.join(outdir, f"filtered_{file}")
        with open(output_file, 'w') as out_handle:
            for species, stop_dict in all_orfs.items():
                for stop_codon, orf_list in stop_dict.items():
                    # 按长度降序排序，取第一个（最长）
                    orf_list.sort(key=lambda x: x[0], reverse=True)
                    longest_orf = orf_list[0][1]
                    SeqIO.write(longest_orf, out_handle, 'fasta')

关键修正点： #

1. 正确获取标识符：用record.id.strip()去除FASTA标题前的空格，避免分割错误。
2. 长度转为整数：确保比较大小的正确性。
3. 简化存储结构：用列表存储每个组的ORF（长度+记录），后续直接排序取最大值，比嵌套字典更高效。
4. 逻辑顺序调整：先遍历所有记录完成数据存储，再统一筛选最长ORF，避免重复计算。
5. 输出结果：将筛选后的最长ORF写入输出目录下的对应文件，方便后续使用。

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