如何通过区域列表文件提取FASTA文件中用户指定的序列区域

前提文件示例 #

1. test.fasta 内容 #

>seq1
ATCGATCGATCGATCGATCG
>seq2
GGCCGGCCGGCCGGCCGGCC
>seq3
TTAAGGCCGGCCGGTTAACC

2. range.txt 内容（格式：序列ID 起始位置 终止位置，FASTA序列位置默认从1开始计数） #

seq1 3 10
seq2 5 15
seq3 2 8

预期输出 #

>seq1:3-10
CGATCGAT
>seq2:5-15
CCGGCCGGCCG
>seq3:2-8
TAAGGCC

可靠实现方案 #

以下提供两种适配Bash环境的解决方案，覆盖无额外工具依赖和专业生物信息工具两种场景：

方案1：纯awk实现（无需额外安装工具） #

该脚本支持处理多行拆分的FASTA序列，严格遵循1-based位置规则：

awk '
BEGIN {
    # 预读取range.txt的提取规则，存入数组
    while ((getline < "range.txt") > 0) {
        ranges[$1] = $2 "," $3
    }
    close("range.txt")
    current_id = ""
    seq = ""
}
/^>/ {
    # 处理已缓存的上一条序列（如果存在未完成的提取）
    if (need_extract && length(seq) >= end) {
        print substr(seq, start, end - start + 1)
    }
    # 重置当前状态
    current_id = substr($0, 2)
    seq = ""
    need_extract = 0
    # 检查当前ID是否在提取规则中
    if (current_id in ranges) {
        split(ranges[current_id], pos, ",")
        start = pos[1]
        end = pos[2]
        need_extract = 1
        print ">" current_id ":" start "-" end
    }
    next
}
need_extract {
    # 拼接多行序列
    seq = seq $0
    # 序列长度足够时立即提取并重置
    if (length(seq) >= end) {
        print substr(seq, start, end - start + 1)
        seq = ""
        need_extract = 0
    }
}
END {
    # 处理文件末尾的最后一条序列
    if (need_extract && length(seq) >= end) {
        print substr(seq, start, end - start + 1)
    }
}
' test.fasta

方案2：使用seqtk（生物信息学专用工具，需提前安装） #

seqtk是处理FASTA/Q文件的轻量高效工具，命令更简洁：

# 将range.txt转换为seqtk支持的格式（每行：ID:起始-终止）
awk '{print $1 ":" $2 "-" $3}' range.txt > extract_rules.txt
# 执行序列提取
seqtk subseq test.fasta extract_rules.txt

注：seqtk可通过包管理器安装，例如apt install seqtk（Debian/Ubuntu）或brew install seqtk（macOS）。

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