使用Awk按条件过滤文件时的范围匹配问题求助

输入文件说明 #

1. gene_positions.txt #

共3列，分别为Chromosome（染色体）、Gene_start（基因起始位置）、Gene_end（基因终止位置），示例（仅4号染色体内容）：

Chromosome Gene_start Gene_end 
4 25121627 25167204
4 170981373 171017850
4 37592422 37692998
4 84382094 84411290

2. gwas_results.txt #

共14列，核心关注前两列：CHROM（染色体编号）、POS（位点位置），示例：

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中满足以下条件的行并输出到lookup_output.txt：

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

错误脚本及问题 #

使用的Awk脚本如下：

#!/bin/bash
#SBATCH --job-name=genelookup
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --cpus-per-task=1
#SBATCH --mem=16gb
#SBATCH --time=01:00:00
#SBATCH --export=NONE

cd /data/genome/hj/GWAS_results 

awk '
NR == FNR {
   start[$1] = $2
   end[$1] = $3
   next
}
$1 in start && $2 >= start[$1] && $3 <= end[$1]
' gene_positions.txt gwas_results.txt > lookup_output.txt

问题：输出结果包含不符合区间的行，例如4号染色体的POS=84529777，超出了gene_positions.txt中4号染色体的最大区间（84382094-84411290）。原因是原脚本用染色体号作为数组键，会覆盖同一染色体的多个区间，只保留最后一行的start和end值，导致无法匹配所有区间。

正确解决方案 #

修改Awk脚本，为每个染色体存储所有区间而非仅最后一个，脚本如下：

#!/bin/bash
#SBATCH --job-name=genelookup
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --cpus-per-task=1
#SBATCH --mem=16gb
#SBATCH --time=01:00:00
#SBATCH --export=NONE

cd /data/genome/hj/GWAS_results 

awk '
NR == FNR {
    # 跳过gene_positions.txt的表头
    if (NR == 1) next
    # 为每个染色体存储所有区间，格式为"start1,end1;start2,end2;"
    chrom = $1
    intervals[chrom] = intervals[chrom] $2 "," $3 ";"
    next
}
# 处理gwas_results.txt
{
    chrom = $1
    pos = $2
    # 保留并输出gwas_results.txt的表头
    if (chrom == "CHROM") {
        print $0
        next
    }
    # 当前染色体无匹配区间则跳过
    if (!(chrom in intervals)) next
    # 遍历该染色体的所有区间
    n = split(intervals[chrom], arr, ";")
    for (i=1; i<n; i++) {
        split(arr[i], se, ",")
        start = se[1]
        end = se[2]
        if (pos >= start && pos <= end) {
            print $0
            break  # 匹配到一个区间即输出，避免重复输出同一行
        }
    }
}
' gene_positions.txt gwas_results.txt > lookup_output.txt

脚本说明 #

1. 处理gene_positions.txt时，为每个染色体存储所有区间，用分号分隔不同区间，逗号分隔每个区间的起始和终止位置
2. 处理gwas_results.txt时，保留并输出表头；对每个位点，遍历对应染色体的所有区间，判断POS是否落在任意区间内，满足条件则输出该行
3. 若不需要保留GWAS结果的表头，可删除对应表头处理的代码块

预期输出示例 #

仅输出CHROM匹配且POS落在对应Gene_start和Gene_end区间内的行：

4   25121630    4:84392377:G:T  G   T   T   ADD 1635    -0.00336762 0.0415539   -0.0848117  0.0780765   -0.0810422  0.935418
4   25131627    4:84392378:A:T  A   T   T   ADD 1635    -0.00336762 0.0415539   -0.0848117  0.0780765   -0.0810422  0.935418
4   170981374   4:84415536:CA:C CA  C   C   ADD 1635    -0.0359081  0.0576552   -0.14891    0.0770941   -0.622807   0.533499
4   171017750   4:84415882:AAACATG:A    AAACATG A   A   ADD 1635    -0.00552875 0.0522986   -0.108032   0.0969745   -0.105715   0.915821
4   37492998    4:84454092:T:TTATATATG  T   TTATATATG   TTATATATG   ADD 1635    0.0148902   0.0442058   -0.0717516  0.101532    0.336838    0.736283
4   37592998    4:84458776:G:GTA    G   GTA GTA ADD 1635    -0.000157512    0.0252463   

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