生物信息学新手求助：如何读取基因对列表并生成对应FASTA文件

核心思路 #

先将大型FASTA文件中的基因ID与对应序列存入字典（实现快速查找），再逐行读取基因对列表，为每一对基因生成独立的FASTA文件。

完整Python脚本 #

# 1. 读取FASTA文件，构建基因-序列映射字典
gene_seq_dict = {}
current_gene = None
current_seq = []

with open('genes.faa', 'r') as faa_file:
    for line in faa_file:
        line = line.strip()
        if not line:
            continue
        # 识别基因ID行
        if line.startswith('>'):
            # 保存上一个基因的序列（如果存在）
            if current_gene is not None:
                gene_seq_dict[current_gene] = ''.join(current_seq)
            # 更新当前基因ID，重置序列缓存
            current_gene = line.lstrip('>')
            current_seq = []
        # 收集序列片段（兼容单行/多行序列）
        else:
            current_seq.append(line)
# 保存最后一个基因的序列
if current_gene is not None:
    gene_seq_dict[current_gene] = ''.join(current_seq)

# 2. 处理基因对列表，生成对应FASTA文件
pair_index = 1
with open('gene.pairs.txt', 'r') as pairs_file:
    for line in pairs_file:
        line = line.strip()
        if not line:
            continue
        # 分割制表符分隔的基因对
        gene1, gene2 = line.split('\t')
        # 生成输出文件名
        output_file = f'pair_{pair_index}.fasta'
        # 写入FASTA内容
        with open(output_file, 'w') as out:
            out.write(f'>{gene1}\n{gene_seq_dict[gene1]}\n')
            out.write(f'>{gene2}\n{gene_seq_dict[gene2]}\n')
        print(f"已生成文件: {output_file}")
        pair_index += 1

脚本细节说明 #

• FASTA读取逻辑：兼容单行或多行序列的真实FASTA格式，用列表缓存序列片段后拼接，避免因序列换行导致的错误。
• 基因对处理逻辑：逐行读取基因对文件，通过字典快速调取序列，无需反复读取大型FASTA文件，提升运行效率。
• 容错扩展：若担心基因ID不存在，可在写入前添加判断（跳过缺失基因的对或打印警告）：

  if gene1 not in gene_seq_dict or gene2 not in gene_seq_dict:
      print(f"警告：基因对{gene1}/{gene2}存在未找到的基因，跳过")
      continue
  

使用步骤 #

1. 将脚本保存为generate_pair_fastas.py
2. 确保genes.faa、gene.pairs.txt与脚本在同一目录
3. 运行脚本：python generate_pair_fastas.py

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