如何从FASTA文件中随机删除指定比例核苷酸并修复脚本报错

原有shell脚本存在两层核心问题，无法实现需求：

• 逻辑适配错误：参考的是文件按行随机删除方案，但FASTA格式中序列可能拆分为多行存储，即便序列为单行，sed删除行的操作是移除整行内容，无法精准删除单个核苷酸字符，完全不匹配核苷酸级别的删除需求
• 语法报错原因：shuf生成的待删除位置是换行分隔的数字列表，直接传入printf作为参数时，zsh会尝试将未正确包裹的数字串解析为算术表达式，缺少运算符就会抛出bad math expression报错。即便绕过该报错，sed按行删除的逻辑也无法得到正确结果。

处理百万级长度的核苷酸序列，用Python配合Biopython实现是最稳定高效的选择，运行耗时不到1秒，自动兼容所有标准FASTA格式，不会出现格式解析错误。

首先安装依赖：
pip install biopython

然后使用以下脚本，修改开头的文件路径配置即可运行：

import random
from Bio import SeqIO
from Bio.Seq import Seq

# 配置参数
input_fasta = "/Users/home/DETECTION/GCA_900186885.1_48903_D01_genomic_reformatted.fa"
output_configs = [
    {"keep_ratio": 0.9, "output_path": "delete_10percent.fasta"},
    {"keep_ratio": 0.85, "output_path": "delete_15percent.fasta"},
    {"keep_ratio": 0.8, "output_path": "delete_20percent.fasta"}
]
random_seed = 42  # 固定随机种子，结果可复现，不需要可注释该行
line_width = 60  # 输出FASTA每行核苷酸长度，符合通用格式规范

random.seed(random_seed)
# 2M长度序列直接读入内存无压力，超大文件可改为逐序列流式处理
records = list(SeqIO.parse(input_fasta, "fasta"))

for config in output_configs:
    keep_ratio = config["keep_ratio"]
    out_path = config["output_path"]
    processed_records = []
    for rec in records:
        seq = str(rec.seq)
        seq_len = len(seq)
        keep_num = round(seq_len * keep_ratio)
        # 生成保留位置索引并排序，保证核苷酸原有顺序不变
        keep_pos = sorted(random.sample(range(seq_len), keep_num))
        new_seq = "".join([seq[i] for i in keep_pos])
        rec.seq = Seq(new_seq)
        processed_records.append(rec)
    # 按格式写入输出文件
    total_len = 0
    with open(out_path, "w") as f:
        for rec in processed_records:
            f.write(f">{rec.id} {rec.description}\n")
            seq_str = str(rec.seq)
            total_len += len(seq_str)
            for i in range(0, len(seq_str), line_width):
                f.write(seq_str[i:i+line_width] + "\n")
    print(f"已生成{out_path}，总核苷酸长度：{total_len}")

使用说明 #

• 脚本默认固定随机种子为42，每次运行生成的删除结果一致，方便实验复现；如果需要每次生成不同的随机结果，注释掉random.seed(random_seed)行即可
• 默认输出FASTA按每行60个核苷酸折行，符合通用FASTA格式规范，可通过修改line_width参数调整折行长度
• 脚本运行后会打印每个输出文件的实际核苷酸总长度，可直接和预期值核对。如果需要严格匹配你给出的固定长度（1918928/1812321/1705714），可以把对应配置下计算keep_num的行替换为固定值赋值，比如10%删除比例的配置中直接写keep_num = 1918928即可
• 不建议用纯shell工具实现该需求：如果把每个核苷酸拆成单行再做删除，200多万行的文本处理性能极差，还容易因为参数长度限制、FASTA头误处理等问题得到错误结果。

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