如何用Bash提取子集基因FASTA序列的2-7位并格式化输出

阿华AIGC实验室

2026-5-9

提取FASTA序列指定片段并调整输出格式的Bash方案

你已经用grep -w -A 1 -f subset.txt all_genes.fasta拿到了目标基因的FASTA条目，接下来只需要用awk来处理这些输出，就能轻松提取序列的2-7位并调整格式。

基础解决方案（序列为单行的情况）

直接把grep的输出管道给awk处理，一行命令搞定：

grep -w -A 1 -f subset.txt all_genes.fasta | awk '
/^>/ { id = substr($0, 2); next }  # 捕获基因ID，去掉开头的>
{ print substr($0, 2, 6) "\t" id } # 提取序列第2-7位，和ID按要求格式输出
'

命令解释：

grep部分：你已经熟悉的逻辑——匹配子集文件中的基因ID，输出对应的ID行（>开头）和紧随其后的序列行
awk部分：
- 当遇到以>开头的行时，用substr($0,2)去掉开头的>，把剩下的基因ID存到变量id里，然后跳到下一行处理
- 处理序列行时，substr($0,2,6)表示从第2个字符开始，截取6个字符（正好对应2-7位），然后用制表符分隔，跟上之前存的基因ID输出

用你的示例数据测试，会得到完全符合期望的输出：

GAGGUA	mmu-let-7g-5p MIMAT0000121
GAGGUA	mmu-let-7i-5p MIMAT0000122

进阶：处理多行序列的情况

如果你的FASTA文件里存在序列换行的情况（比如长序列被拆成了多行），上面的命令就会失效，这时候需要先把每个基因的序列合并成单行，再提取片段。可以用下面的awk命令直接处理原始FASTA文件，不需要先调用grep：

awk '
NR == FNR { ids[$0] = 1; next }  # 读取子集文件，把基因ID存到数组里
/^>/ {
    if (current_id != "" && ids[current_id]) {
        # 输出上一个匹配基因的结果
        print substr(seq, 2, 6) "\t" current_id
    }
    current_id = substr($0, 2)
    seq = ""
    next
}
ids[current_id] { seq = seq $0 }  # 合并当前基因的所有序列行
END {
    # 处理最后一个匹配的基因
    if (current_id != "" && ids[current_id]) {
        print substr(seq, 2, 6) "\t" current_id
    }
}
' subset.txt all_genes.fasta

这个命令会先读取子集文件的所有ID，然后遍历FASTA文件，把每个匹配基因的序列合并成单行，最后提取2-7位并输出。