如何用awk实现双文件匹配 筛选染色体一致的基因区间CpG位点

核心问题说明 #

你原有命令的错误根源是存储位点时仅用位置作为数组键，没有关联对应的染色体号，无法做染色体一致性校验，且200万行位点的情况下遍历逻辑效率极低，很容易出现性能瓶颈。

最优修改方案（推荐，性能更高） #

优先读取文件更小的基因注释文件，按染色体分组存储区间，处理位点时仅校验同染色体的区间，时间复杂度大幅降低：

awk '
# 读取第一个输入文件：基因注释文件
FNR == NR {
    # 跳过注释文件的表头
    if (FNR == 1) next
    chr = $2
    # 按染色体分组存储基因的起始、终止位置和完整信息
    gene_start[chr][FNR] = $3
    gene_end[chr][FNR] = $4
    gene_detail[chr][FNR] = $0
    next
}
# 读取第二个输入文件：CpG位点文件
{
    current_chr = $1
    current_pos = $2
    # 仅当前染色体有注释时才做匹配，跳过无注释的染色体
    if (current_chr in gene_start) {
        # 遍历当前染色体下的所有基因区间
        for (idx in gene_start[current_chr]) {
            # 可根据需求调整为开区间判断 > / <
            if (current_pos >= gene_start[current_chr][idx] && current_pos <= gene_end[current_chr][idx]) {
                # 输出格式：位点完整信息 + 对应基因完整信息，可自行调整输出字段
                print $0, gene_detail[current_chr][idx]
            }
        }
    }
}' 基因注释文件路径 CpG位点文件路径

调用时请把命令里的文件路径替换成你实际的文件路径即可。

兼容原有文件顺序的修改方案 #

如果你一定要保持原有“位点文件在前、注释文件在后”的输入顺序，可以用以下修改版，不过性能会差一些：

awk '
# 读取第一个输入文件：CpG位点文件
FNR == NR {
    # 用“染色体号:位点位置”作为数组键，关联存储位点完整信息
    site_key = $1 ":" $2
    site_info[site_key] = $0
    next
}
# 读取第二个输入文件：基因注释文件，跳过表头
FNR > 1 {
    curr_chr = $2
    curr_start = $3
    curr_end = $4
    # 遍历所有位点，先匹配染色体再判断区间
    for (key in site_info) {
        split(key, key_arr, ":")
        if (key_arr[1] == curr_chr && key_arr[2] > curr_start && key_arr[2] < curr_end) {
            print site_info[key], $0
        }
    }
}' file1.txt file.txt

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