如何使用awk匹配染色体区间 为测序覆盖度文件新增基因名列

原代码问题说明 #

• 未使用FNR==NR区分双文件处理阶段，无法提前将基因区间数据存入内存备用
• awk 不存在file1$1/file2$1这类直接跨文件引用字段的语法，必须先将第一个文件的内容存入自定义数组，处理第二个文件时才能调用

核心实现逻辑 #

采用你提到的FNR==NR双文件处理语法：先读取存储基因区间的File2，将所有基因的染色体、起始位置、终止位置、基因名存入数组；再遍历测序覆盖度文件File1，逐行判断位点是否落在已存储的基因区间内，匹配成功则在原行后追加对应的基因名。

完整可执行代码 #

awk '
# 处理第一个输入文件：基因区间文件File2
FNR == NR {
    chr = $1
    # 按染色体分组存储基因的终止位置和名称
    gene_end[chr][$2] = $3
    gene_name[chr][$2] = $4
    next
}
# 处理第二个输入文件：覆盖度文件File1
{
    cur_chr = $1
    pos = $2
    matched = ""
    # 遍历当前染色体下所有基因区间做匹配
    for (start in gene_end[cur_chr]) {
        if (pos >= start && pos <= gene_end[cur_chr][start]) {
            matched = gene_name[cur_chr][start]
            break
        }
    }
    # 输出原内容+匹配到的基因名，未匹配则追加空值
    print $0, matched
}' File2 File1 > 结果文件.tsv

注意事项 #

• 输入文件的顺序不能颠倒，必须先传基因区间文件File2，再传覆盖度文件File1
• 若需要保持原文件的制表符分隔格式，可将最后一行输出语句修改为print $0 "\t" matched
• 若需要过滤掉未匹配到基因的位点，可在输出前添加判断条件：if (matched != "") print $0, matched

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