Python实现随机向FASTA基因组插入指定序列并记录位置

错误原因分析 #

最初出现的TypeError: can only concatenate str (not "list") to str，核心原因是：

• hundred_seqs是SeqRecord对象的列表，并非字符串类型，直接传入insert函数会尝试将基因组字符串与列表拼接，触发类型不匹配错误
• 即便将序列转为字符串，直接传入列表也不符合“插入到随机位置”的需求，正确逻辑应为逐个序列插入到独立位置

核心需求实现方案 #

要完成「筛选指定长度序列、无重叠随机插入、记录插入位置」三个核心目标，关键步骤如下：

1. 筛选400-500bp的FASTA序列，用random.sample随机选取100个（避免重复选取同一序列）
2. 生成100个不重叠的随机插入位置：通过降序排序位置实现无重叠插入（从后往前插入，避免前面的插入操作影响后续位置的索引）
3. 插入过程中记录每个序列的关键信息（ID、长度）及插入位置（原始基因组中的位置、最终基因组中的起始位置）

完整可运行代码 #

from Bio import SeqIO
import random

# 1. 导入并筛选目标序列（400-500bp）
def filter_sequences(input_fasta):
    # 读取FASTA文件为SeqRecord对象列表
    all_seqs = list(SeqIO.parse(input_fasta, "fasta"))
    # 筛选长度在400-500之间的序列
    filtered_seqs = [seq for seq in all_seqs if 400 < len(seq) < 500]
    # 随机选取100个（无重复）
    return random.sample(filtered_seqs, k=100)

# 2. 读取目标基因组FASTA
def load_genome(genome_fasta):
    # 读取第一条序列（假设基因组为单条序列）
    genome_record = next(SeqIO.parse(genome_fasta, "fasta"))
    return genome_record.id, str(genome_record.seq)

# 3. 无重叠插入序列并记录位置
def insert_sequences(genome_seq, insert_seqs):
    genome = list(genome_seq)  # 转为列表提升插入效率
    insert_records = []  # 存储插入记录：[(序列ID, 原始位置, 最终起始位置, 序列长度), ...]
    original_length = len(genome)

    # 生成100个随机原始插入位置（范围0到原始基因组长度）
    insert_positions = [random.randint(0, original_length) for _ in range(len(insert_seqs))]
    # 将位置与序列配对，按位置降序排序（从后往前插入，避免位置偏移）
    sorted_pairs = sorted(zip(insert_positions, insert_seqs), key=lambda x: x[0], reverse=True)

    for original_pos, seq in sorted_pairs:
        seq_str = str(seq.seq)
        seq_len = len(seq_str)
        # 计算最终起始位置：原始位置 + 之前插入的所有序列总长度
        inserted_total = sum(record[3] for record in insert_records)
        final_start_pos = original_pos + inserted_total
        # 插入序列到基因组
        genome[final_start_pos:final_start_pos] = list(seq_str)
        # 记录信息
        insert_records.append((seq.id, original_pos, final_start_pos, seq_len))

    # 将列表转回字符串
    new_genome = ''.join(genome)
    # 反转记录（恢复原始随机顺序）
    insert_records.reverse()
    return new_genome, insert_records

# 4. 保存结果
def save_results(genome_id, new_genome, insert_records, output_genome_path, output_log_path):
    # 保存修改后的基因组FASTA（按标准格式换行）
    with open(output_genome_path, "w") as f:
        f.write(f">{genome_id}\n")
        for i in range(0, len(new_genome), 60):
            f.write(new_genome[i:i+60] + "\n")
    # 保存插入位置日志
    with open(output_log_path, "w") as f:
        f.write("序列ID\t原始插入位置\t最终起始位置\t序列长度\n")
        for record in insert_records:
            f.write(f"{record[0]}\t{record[1]}\t{record[2]}\t{record[3]}\n")

# 主执行流程
if __name__ == "__main__":
    # 配置文件路径
    target_seqs_path = "hs_sequences.txt"
    genome_path = "hg37_chr1-1000l.fna"
    output_genome = "modified_genome.fasta"
    output_log = "insertion_log.txt"

    # 执行操作
    selected_seqs = filter_sequences(target_seqs_path)
    genome_id, genome_seq = load_genome(genome_path)
    new_genome, insert_records = insert_sequences(genome_seq, selected_seqs)
    save_results(genome_id, new_genome, insert_records, output_genome, output_log)

    print("插入完成！已生成修改后的基因组文件和插入位置日志")

代码关键说明 #

• 序列筛选：用random.sample替代循环random.choice，避免重复选取同一序列
• 无重叠插入：将插入位置按降序排序后从后往前插入，确保每次插入的位置不受之前操作影响，彻底避免序列重叠
• 位置记录：完整存储每个序列的ID、原始插入位置、最终起始位置及长度，方便后续验证与分析
• 高效操作：将基因组转为列表进行插入，相比字符串拼接大幅提升效率（字符串为不可变对象，频繁拼接会产生大量临时对象）

原代码问题修正 #

• 修复了原代码中未定义变量（如s、retro）的问题
• 替换了原代码中“用逗号拼接所有序列插入到单一位置”的错误逻辑，实现了每个序列插入到独立随机位置的需求

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