多行序列文本中GCAGC/GCTGC匹配位置查找及流式实现方法

该方案无需加载整条序列到内存，仅保留最多4个碱基的跨行缓存，内存占用恒定，可支持任意大小的输入文件。

处理脚本 #

awk '
/^>/ {
    if (seq_name != "") {
        prev_suffix = ""
        offset = 0
    }
    seq_name = substr($0, 2)
    next
}
{
    current_seq = prev_suffix $0
    for (i=1; i <= length(current_seq) - 4; i++) {
        kmer = substr(current_seq, i, 5)
        if (kmer == "GCAGC" || kmer == "GCTGC") {
            actual_pos = offset + i - length(prev_suffix)
            print seq_name, actual_pos, kmer
        }
    }
    offset += length($0)
    prev_suffix = substr($0, length($0) - 3)
}
' 你的输入序列文件路径

逻辑说明 #

• 遇到>开头的行时自动提取序列ID，同时重置缓存和偏移计数，天然支持多序列文件处理
• 每次处理序列行时，会把上一行末尾的4个碱基拼接到当前行开头，避免跨行的匹配位点被遗漏
• 遍历每个可能的5碱基片段，匹配目标模式后计算全局1-based起始位置，符合生物学序列计数习惯
• 运行全程仅保留当前序列ID、全局偏移量、4个碱基的后缀缓存三个变量，内存占用不随文件大小增长

输出说明 #

输出完全符合要求的三列格式：第一列为序列ID，第二列为匹配起始位置，第三列为匹配到的具体模式，重叠匹配位点会单独输出。

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