使用AWK命令匹配过滤两文件数据时无输出的问题排查

问题描述 #

现有两个文件：

1. gene_positions.txt（3列：染色体、基因起始位置、基因结束位置）：

>> cat gene_positions.txt
Chromosome Gene_start Gene_end 
15 26788693 27189686
1 22004791 22115099
20 44689624 44723591
20 61872136 61909046

2. gwas_results.txt（14列，仅关注前两列CHROM、POS）：

>> cat gwas_results.txt
CHROM   POS ID  REF ALT A1  TEST    OBS_CT  BETA    SE  L95 U95 T_STAT  P
1   22004791    1:54712:T:TTTTC T   TTTTC   T   ADD 1597    -0.00256645 0.0313666   -0.0640439  0.058911    -0.0818212  0.934799
1   22105000    1:54712:T:TC    T   TTTTC   T   ADD 1597    -0.00259875 0.0313666   -0.0640439  0.058911    -0.0818216  0.94799
1   825410  rs13303179:825410:G:A   G   A   G   ADD 1597    -0.017462   0.0235123   -0.0635454  0.0286213   -0.742676   0.457788

需求：筛选出gwas_results.txt中满足以下条件的行并输出到test_output_results.txt：

• CHROM列与gene_positions.txt的Chromosome列精确匹配
• POS列数值落在对应染色体的任意一个Gene_start与Gene_end区间内

使用以下AWK命令后得到0输出：

awk 'NR==FNR {chromosome[$1]; start[$1]; end[$1]; next} 
   ($1 in chromosome) && ($2 >= start[$1]) && ($2 <= end[$1]) {print}'\
 gene_positions.txt gwas_results.txt > test_output_results.txt

预期输出：

1   22004791    1:54712:T:TTTTC T   TTTTC   T   ADD 1597    -0.00256645 0.0313666   -0.0640439  0.058911    -0.0818212  0.934799
1   22105000    1:54712:T:TC    T   TTTTC   T   ADD 1597    -0.00259875 0.0313666   -0.0640439  0.058911    -0.0818216  0.94799

错误原因 #

1. 数组赋值无效：原脚本中chromosome[$1]; start[$1]; end[$1];仅创建了数组键，但未将对应行的Gene_start、Gene_end值赋值给数组，导致start[$1]和end[$1]为空，无法进行有效数值比较。
2. 多区间染色体处理错误：同一染色体可能对应多个基因区间（如示例中20号染色体有两个区间），原脚本用染色体号作为数组键会覆盖之前的区间值，无法匹配所有区间。
3. 未跳过表头：两个文件的第一行都是表头字符串，会被误处理为数据行，干扰染色体匹配和数值比较逻辑。

正确的AWK命令 #

脚本版本（可读性更高） #

创建脚本文件match_intervals.awk：

# 处理基因位置文件，跳过表头，存储每个染色体的所有区间
NR==FNR {
    if (NR == 1) next
    # 用分号分隔同一染色体的多个区间，每个区间用逗号分隔起始/结束位置
    chrom_intervals[$1] = chrom_intervals[$1] $2 "," $3 ";"
    next
}

# 处理GWAS结果文件，跳过表头
NR != FNR {
    if (NR == FNR + 1) next
    chrom = $1
    pos = $2
    # 检查当前染色体是否有对应的区间
    if (chrom in chrom_intervals) {
        # 拆分该染色体的所有区间
        split(chrom_intervals[chrom], intervals, ";")
        # 遍历每个区间进行匹配
        for (i in intervals) {
            if (intervals[i] == "") continue
            split(intervals[i], se, ",")
            start_pos = se[1]
            end_pos = se[2]
            # 匹配成功则输出该行并跳出循环（避免重复输出）
            if (pos >= start_pos && pos <= end_pos) {
                print
                break
            }
        }
    }
}

运行命令：

awk -f match_intervals.awk gene_positions.txt gwas_results.txt > test_output_results.txt

单行版本（直接执行） #

awk 'NR==FNR{if(NR==1)next;chrom_intervals[$1]=chrom_intervals[$1]$2","$3";";next}NR!=FNR{if(NR==FNR+1)next;chrom=$1;pos=$2;if(chrom in chrom_intervals){split(chrom_intervals[chrom],intervals,";");for(i in intervals){if(intervals[i]=="")continue;split(intervals[i],se,",");start=se[1];end=se[2];if(pos>=start&&pos<=end){print;break}}}}' gene_positions.txt gwas_results.txt > test_output_results.txt

结果验证 #

执行命令后，test_output_results.txt会生成预期的两行输出，完全符合需求。

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