如何编写Python脚本计算MySequences.fasta中序列的反向互补序列？

需求：编写可自定义命名的Python脚本，分析FASTA文件MySequences.fasta，输出其中每条序列的反向互补序列。

原代码存在的问题 #

原代码有几个关键问题导致无法正常工作：

• 第一次遍历文件后，文件指针已移至末尾，第二次遍历不会执行任何操作
• 仅保存了最后一行序列，未处理FASTA序列跨多行的情况
• 导入了itertools.repeat但未实际使用
• 未关联输出反向互补序列对应的标题

修正后的完整代码 #

你可以将脚本保存为任意名称（比如fasta_rc_generator.py），代码如下：

def reverse_complement(seq):
    # 用字典实现碱基互补映射，比多条件判断更高效
    complement_map = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'}
    # 生成互补序列后反转，得到反向互补
    return ''.join([complement_map[base] for base in seq[::-1]])

def process_fasta(filename):
    sequences = {}
    current_header = None
    current_seq = []
    
    with open(filename, 'r') as file:
        for line in file:
            line = line.strip()
            if not line:
                continue
            # 识别标题行
            if line.startswith('>'):
                # 如果已有未保存的序列，先存入字典
                if current_header and current_seq:
                    sequences[current_header] = ''.join(current_seq)
                current_header = line
                current_seq = []
            else:
                # 收集序列行（处理多行序列）
                current_seq.append(line)
        # 处理最后一条序列
        if current_header and current_seq:
            sequences[current_header] = ''.join(current_seq)
    
    # 输出每条序列的反向互补
    for header, seq in sequences.items():
        rc_seq = reverse_complement(seq)
        print(f"{header}_reverse_complement")
        print(rc_seq)
        print()  # 空行分隔不同序列

if __name__ == "__main__":
    # 可自定义文件名，这里默认处理MySequences.fasta
    input_file = "MySequences.fasta"
    process_fasta(input_file)

代码说明 #

• 反向互补函数：使用字典映射碱基互补关系，先反转序列再生成互补，逻辑更清晰高效
• FASTA文件处理：使用with语句自动管理文件打开/关闭，正确收集跨多行的序列，并将标题与序列一一对应存入字典
• 自定义设置：你可以修改input_file变量指定目标FASTA文件，也可以将输出结果写入新文件（替换print部分为文件写入代码即可）
• 可扩展性：如果需要处理大小写混合的碱基（比如小写a/t/c/g），可以在映射字典中添加对应键值，或者先将序列转为大写

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