如何编写Bash脚本处理FASTA文件：修改标题行并去除N碱基

需要编写Bash脚本处理基因组FASTA文件，完成以下操作：

• 对所有以>开头的标题行，提取"chromosome"后的数字，替换为>chr{digit}格式（如>chr1），删除该行其他内容
• 对非标题行（序列行），移除所有"N"（包括大小写），并将同一条染色体的序列拼接为单行
• 原尝试的两段代码运行超1小时仍未完成，输入输出示例如下：

输入示例：

>NC_000001.11 Homo sapiens chromosome 1, GRCh38.p14 Primary Assembly
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

gaaagtGAGGTTGCCTGCCCTGTCTCCTA
>NC_000001.12 Homo sapiens chromosome 2, GRCh38.p14 Primary Assembly
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
gaaagtGAGGTTGCCTGCCCTGTCTCCTAGAGGTTGCCTGCCCTGTCTCCTA

期望输出：

>chr1 
gaaagtGAGGTTGCCTGCCCTGTCTCCTA

>chr2
gaaagtGAGGTTGCCTGCCCTGTCTCCTAGAGGTTGCCTGCCCTGTCTCCTA

原尝试代码：
处理标题行（存在语法错误，无输入文件）：

grep '^>' | grep -oP "chromosome \s+\K\w+" | cat "chr" GCF_000001405.40_GRCh38.p14_genomic.fna

处理序列行（多管道导致效率低下）：

awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' < GCF_000001405.40_GRCh38.p14_genomic.fna | sed 's/N//g' | tr "\t" "\n" > genome.fna

用单段awk脚本一次性完成所有操作，避免重复读取文件和多管道开销，大幅提升处理速度：

脚本代码（保存为process_genome.awk） #

/^>/ {
    # 输出格式化后的标题行
    if (match($0, /chromosome[[:space:]]+([0-9]+)/, arr)) {
        print ">chr" arr[1]
    }
    # 输出上一条染色体的有效序列（如果存在）
    if (seq != "") {
        print seq "\n"
        seq = ""
    }
    next
}
{
    # 移除当前行的所有N/n，拼接有效序列
    gsub(/[Nn]/, "")
    seq = seq $0
}
END {
    # 输出最后一条染色体的有效序列
    if (seq != "") print seq
}

执行命令 #

awk -f process_genome.awk GCF_000001405.40_GRCh38.p14_genomic.fna > genome.fna

关键说明 #

1. 标题行处理：通过正则精准匹配chromosome后的数字，直接输出>chr{digit}格式的标题，同时输出上一条染色体的拼接序列。
2. 序列处理：用gsub移除所有大小写N，将同一条染色体的所有有效序列片段拼接为单行，自动忽略全N的空行。
3. 效率优化：仅读取一次输入文件，无额外管道进程，内存占用低，处理大FASTA文件速度远超原方案。

原代码的问题 #

1. 第一段grep未指定输入文件，会一直等待标准输入，导致脚本卡住。
2. 分两次处理文件，重复读取大基因组FASTA，开销翻倍。
3. 多管道（awk→sed→tr）会启动多个进程并传递海量数据，大幅降低处理效率。

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