如何通过命令行匹配覆盖指定基因坐标范围的分块vcf.gz文件

完全可以通过Bash命令行工具批量实现匹配，效率远高于人工操作，以下是可直接复用的实现方案：

先把所有vcf分块文件的文件名解析成「染色体、分块起始坐标、分块终止坐标、完整文件名」的结构化索引，再遍历基因坐标表的每一行，匹配同染色体下满足「分块起始≤基因起始 且 分块终止≥基因终止」的文件即可。

1. 提取所有vcf分块的坐标信息生成索引文件，执行以下命令：

ls *.vcf.gz | awk -F'_' '{
    # 适配示例文件名格式：main_programme_ver2_chr[A]_[B]_[C]_VEPannot.vcf.gz
    chr = $4
    block_start = $5
    block_end = $6
    print chr, block_start, block_end, $0
}' > vcf_block_index.txt

若文件名前缀的下划线数量有变化，调整$4/$5/$6的字段序号即可。

2. 遍历基因坐标表批量匹配，假设基因坐标文件名为gene_coords.txt，执行以下命令输出匹配结果：

awk 'NR==FNR{
    # 读入vcf分块索引按染色体分组存储
    chr = $1
    start[chr][NR] = $2
    end[chr][NR] = $3
    file[chr][NR] = $4
    next
}
{
    # 匹配同染色体下覆盖基因全范围的分块
    g_chr = $1
    g_start = $2
    g_end = $3
    g_name = $4
    for(i in start[g_chr]){
        if(start[g_chr][i] <= g_start && end[g_chr][i] >= g_end){
            print g_name, file[g_chr][i]
            break
        }
    }
}' vcf_block_index.txt gene_coords.txt > gene2vcf_mapping.txt

最终输出的gene2vcf_mapping.txt第一列为基因名，第二列为对应的vcf分块文件名，可直接用于下游分析。

• 若存在基因跨多个vcf分块的场景，删除上述代码中的break即可输出所有覆盖该基因的分块文件
• 脚本依赖系统原生awk即可运行，1000+分块、2万+基因的匹配耗时在秒级
• 坐标数值默认按整数比较，无需额外转换格式

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