Unix下用映射文件重命名FASTA序列名遇sed报错求助

问题原因 #

你之前用sed -f scaffold_mapping1.txt报错，是因为sed -f要求输入的是sed命令脚本，而你的映射文件是原始序列名称列表，其中的>符号会被sed误识别为未知命令，导致报错。

方法一：生成sed脚本后处理 #

步骤1：将映射文件转换为sed可执行的替换脚本 #

运行以下命令，从scaffold_mapping.txt中提取映射关系，生成sed脚本：

sed -n '2,$s/^>\(scaffold_[0-9]*_c[0-9]*\)\(.*\)/s/^> *\1$/>\1\2/p/' scaffold_mapping.txt > rename_script.sed

• 说明：
  • -n：仅输出经过处理的行
  • 2,$：跳过第一行注释，从第二行开始处理
  • 正则捕获组\(scaffold_[0-9]*_c[0-9]*\)提取短序列名（如scaffold_1_c1），\(.*\)提取后续的_cov...部分
  • 最终生成的sed命令会匹配FASTA中以> （或>）开头的短序列名行，替换为完整名称

步骤2：用生成的脚本处理FASTA文件 #

sed -f rename_script.sed assembly.contigs.fasta > output1.fasta

方法二：用awk直接处理（无需中间脚本） #

awk可以一次性完成映射关系读取和FASTA文件处理，更简洁：

awk '
BEGIN {
    # 读取映射文件，构建短名到完整名的映射表
    while ((getline line < "scaffold_mapping.txt") > 0) {
        if (line ~ /^>scaffold/) {
            split(line, parts, /_cov/);
            short = substr(parts[1], 2);
            full = substr(line, 2);
            map[short] = full;
        }
    }
    close("scaffold_mapping.txt");
}
/^>/ {
    sub(/^> */, "", $0);
    print ">" (map[$0] ? map[$0] : $0);
    next;
}
1' assembly.contigs.fasta > output1.fasta

• 说明：
  • BEGIN块：先读取映射文件，提取短名和对应完整名存入数组map
  • 遇到FASTA的序列名行（以>开头），去掉前缀后查找映射表，输出对应的完整名称；无匹配则保留原名称
  • 非序列名行直接输出（1是awk的简写，代表打印当前行）

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