使用awk批量按CHR匹配与BP窗口规则提取多GWAS文件目标行

现有文件结构 #

• 索引文件：共5000行记录，包含3列，分别为CHR（染色体号）、BP（碱基位置）、FILE（对应GWAS文件名），每行对应一个命名格式为gwas*的独立GWAS文件，结构示例：

CHR     BP          FILE
chr1    12345678    gwas1
chr2    87654321    gwas2
...

• GWAS数据文件：本地共存储5000个与索引FILE列文件名一一对应的gwas*文件，以gwas1为例，单文件包含CHR、BP两列位点数据，结构示例：

CHR     BP                  
chr1    12345678    
chr1    12345679        
chr1    12356777    
...

提取规则 #

针对每个gwas*文件，提取满足以下两个条件的行：

1. 行的CHR列值与索引文件中该文件对应的CHR值完全匹配
2. 行的BP值处于索引中对应BP值上下浮动500000的区间范围内

规则示例：索引中gwas1对应CHR为chr1、BP为12345678，则需提取gwas1中CHR为chr1、BP值介于11845678（BP值减500000）到12845678（BP值加500000）之间的所有行。

现有单文件处理代码 #

目前已实现单文件手动提取的awk逻辑，但未对接索引文件，无法批量处理5000个文件，当前使用环境为GNU Awk 4.0.2，单文件代码如下：

export CHR="chr11"
export BP=107459522
export WINDOW=500000

awk -v CHR=$CHR -v BP_pos=$(($BP + $WINDOW)) -v BP_neg=$(($BP - $WINDOW)) 'BEGIN{FS=OFS="\t"}FNR==1 || ($1 == CHR && $2 < BP_pos && $2 > BP_neg )' gwas1 > gwas1_extract

批量处理实现方案 #

直接运行以下awk命令即可完成全量文件批量处理，无需编写shell循环，处理效率更高：

# 配置区间浮动窗口大小，与原有参数保持一致
WINDOW=500000

awk -v win="$WINDOW" '
BEGIN {
    FS = OFS = "\t"
}
# 读取第一个输入文件（索引文件），存储每个gwas文件的匹配规则
NR == FNR {
    if (FNR > 1) {
        file_target_chr[$3] = $1
        file_bp_low[$3] = $2 - win
        file_bp_high[$3] = $2 + win
    }
    next
}
# 处理后续所有gwas数据文件
{
    current_file = FILENAME
    out_path = current_file "_extract"
    # 每个文件的第一行是表头，直接写入输出文件
    if (FNR == 1) {
        print $0 > out_path
        next
    }
    # 按规则匹配行，符合条件则写入对应输出文件
    if ($1 == file_target_chr[current_file] && $2 > file_bp_low[current_file] && $2 < file_bp_high[current_file]) {
        print $0 >> out_path
    }
}
' 你的索引文件名.txt gwas*

使用说明 #

• 将命令中你的索引文件名.txt替换为实际存储的索引文件名称
• 所有提取结果会自动输出到同目录下，命名规则为原gwas文件名_extract，与单文件处理的输出命名一致
• 代码完全保留原有单文件处理的匹配逻辑，自动保留每个gwas文件的原表头，输出结构符合要求
• 单awk进程处理所有文件，比逐文件调用awk的shell循环方案速度快数倍，适合5000个文件的批量处理场景

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