如何使用Python搜索并修改基因序列中的单个核苷酸

嘿，作为编程新手能写出这样的初始代码已经很棒啦！我来帮你把这个脚本完善，实现你想要的「搜索特定片段+修改单个核苷酸」的功能，一步步来：

先修正现有代码的小问题 #

你的代码已经完成了读取序列的第一步，但有个小bug：number变量的统计逻辑不对，直接用len(GenSeqs)就能拿到正确的序列数量啦，替换掉原来的打印语句就行。另外文件里的header标记应该是>，不是转义后的>，我也一起帮你修正了。

完善核心功能：搜索基序并突变 #

我们来补充用户变量和核心的突变逻辑，每个步骤都写了详细注释，方便你理解：

完整修改后的代码 #

#!/usr/bin/python
from collections import defaultdict
import re, sys, random

# - - - - - H E A D E R - - - - - - - - - - - - - - - - -
'''
Objectives:
1. Read in a sequence
2. Find a specific segment of that sequence
3. Change a letter (mutation)
4. Output the sequence with the mutation
'''
# - - - - - U S E R V A R I A B L E S - - - - - - - -
mssg = " Search and Mutate"
genFile = 'P1.txt'
inFile = 'P1.txt'
inFolder = '.'
# 要搜索的目标基序（示例用GTG，你可以改成自己需要的片段）
site = "GTG"
# 在找到的基序中，要突变的位置（从0开始计数，比如第2个核苷酸填1）
mutation_pos_in_site = 1
# 突变后的核苷酸（示例把原来的碱基改成A）
mutated_nucleotide = "A"
outFile = "Project1-Out.txt"
GenSeqs = defaultdict(lambda: "my own unknown" )
#- - - - - - - - - - - - - - -

# 定义突变函数：输入原始序列、目标基序、基序内突变位置、新核苷酸，返回突变后的序列
def mutate_sequence(sequence, target_site, pos_in_site, new_base):
    # 找到序列中所有匹配基序的起始索引
    match_positions = [match.start() for match in re.finditer(target_site, sequence)]
    
    if not match_positions:
        print(f"⚠️ 警告：在序列中未找到目标基序 '{target_site}'")
        return sequence
    
    # 默认修改第一个匹配的基序，想随机选的话可以改成 random.choice(match_positions)
    start_idx = match_positions[0]
    # 计算突变在整个序列中的绝对位置
    mutation_index = start_idx + pos_in_site
    
    # Python字符串不可直接修改，先转成列表再修改单个字符
    sequence_list = list(sequence)
    sequence_list[mutation_index] = new_base
    
    # 转回字符串返回
    return ''.join(sequence_list)

# - - - - - M A I N - - - - - - - - - - - - - - - - - - - -
print("\n\n", mssg, ". . . . ")

# Task 1: 读取文件中的序列
with open(inFile, 'r') as IN1:  # 用with语句自动关闭文件，更安全
    header = ""
    for line in IN1:
        line = line.strip()
        if line.startswith('>'):
            header = line
        elif header and line:  # 确保有header且当前行是有效序列
            GenSeqs[header] = line

# 打印正确的序列数量
print(f"There are {len(GenSeqs)} gene sequences in file {inFile}")

# Task 2-4: 搜索基序、执行突变、保存结果
mutated_sequences = {}
for seq_header, original_seq in GenSeqs.items():
    print(f"正在处理序列: {seq_header}")
    mutated_seq = mutate_sequence(original_seq, site, mutation_pos_in_site, mutated_nucleotide)
    mutated_sequences[seq_header] = mutated_seq

# 写入输出文件
with open(outFile, 'w') as out_file:
    for header, seq in mutated_sequences.items():
        out_file.write(f"{header}\n")
        out_file.write(f"{seq}\n")

print(f"✅ 突变完成！结果已保存到 {outFile}")

关键部分解释 #

1. 文件读取优化：用with语句管理文件，不用手动关闭；修正了header的匹配逻辑，直接识别文件里的>标记。
2. 突变函数逻辑：
   • 用re.finditer找出所有匹配基序的位置，方便你选择修改哪一个（默认第一个，也可以改成随机选）
   • 因为Python字符串是不可变类型，所以先转成列表修改单个字符，再转回字符串
3. 灵活的用户变量：你只需要修改site、mutation_pos_in_site、mutated_nucleotide这三个变量，就能实现不同的搜索和突变需求。

扩展小技巧 #

• 如果想随机突变序列中的任意单个核苷酸（不基于基序），可以把突变函数改成：

  def random_mutate(sequence, new_base=None):
      if not sequence:
          return sequence
      # 随机选一个序列位置
      random_pos = random.randint(0, len(sequence)-1)
      # 若没指定新核苷酸，随机选ATCG中不同于原碱基的一个
      if not new_base:
          bases = ['A','T','C','G']
          bases.remove(sequence[random_pos])
          new_base = random.choice(bases)
      seq_list = list(sequence)
      seq_list[random_pos] = new_base
      return ''.join(seq_list)
  

• 如果想修改所有匹配的基序，只需要在突变函数里遍历match_positions逐个修改就行。

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