高效合并VCF文件并扩展样本元数据的技术问询

我正在处理一个包含默认列和多个样本相位基因型的大型VCF文件。经其他工具处理后，得到一个变异子集文件，其中包含每个样本的额外元数据。我的目标是将第二个文件的元数据扩展到第一个全量文件的所有变异中，为全量文件的每个变异匹配子集里最近位置的变异的元数据。

面向非生物信息学者的补充说明 #

• 每一行代表一个变异，包含两个核心位置信息：染色体（CHROM）和位点坐标（POS）。
• 扩展规则：为全量VCF的每个变异，在子集VCF中找到同染色体上位置最近的变异，将其元数据（FORMAT列的额外字段、对应样本的元数据值）复制到输出。
• 原全量VCF、子集VCF、输出VCF均为按染色体+坐标排序的文件。

示例文件片段 #

原全量VCF文件： #

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
chr1    1001    rs123   A       G       20.5    PASS    DP=30   GT      0/1
chr1    2001    .       T       C       15.2    PASS    DP=25   GT      1/1
chr1    3001    rs456   C       T       18.7    PASS    DP=28   GT      0/0
chr1    4001    .       G       A       22.0    PASS    DP=35   GT      1/0
chr2    5001    rs789   A       T       25.5    PASS    DP=40   GT      1/1
chr2    6001    .       C       G       30.1    PASS    DP=32   GT      0/1
chr2    7001    rs987   T       A       17.3    PASS    DP=27   GT      1/1
chr3    8001    .       G       C       19.8    PASS    DP=22   GT      0/0
chr3    9001    rs654   A       G       21.4    PASS    DP=26   GT      1/0
chr3    10001   .       T       C       23.2    PASS    DP=29   GT      0/1

处理后子集VCF文件（含额外元数据）： #

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT    Sample1
chr1    1001    rs123   A       G       20.5    PASS    DP=30   GT:X:Y:Z  0/1:0:0:1
chr1    3001    rs456   C       T       18.7    PASS    DP=28   GT:X:Y:Z  0/0:1:3:1
chr2    6001    .       C       G       30.1    PASS    DP=32   GT:X:Y:Z  0/1:2:1:1
chr3    9001    .       T       C       23.2    PASS    DP=29   GT:X:Y:Z  0/1:0:1:0

期望最终输出： #

chr1    1001    rs123   A       G       20.5    PASS    DP=30   GT:X:Y:Z  0/1:0:0:1
chr1    2001    .       T       C       15.2    PASS    DP=25   GT:X:Y:Z  1/1:0:0:1
chr1    3001    rs456   C       T       18.7    PASS    DP=28   GT:X:Y:Z  0/0:1:3:1
chr1    4001    .       G       A       22.0    PASS    DP=35   GT:X:Y:Z  1/0:1:3:1
chr2    5001    rs789   A       T       25.5    PASS    DP=40   GT:X:Y:Z  1/1:2:1:1
chr2    6001    .       C       G       30.1    PASS    DP=32   GT:X:Y:Z  0/1:2:1:1
chr2    7001    rs987   T       A       17.3    PASS    DP=27   GT:X:Y:Z  1/1:2:1:1
chr3    8001    .       G       C       19.8    PASS    DP=22   GT:X:Y:Z  0/0:0:1:0
chr3    9001    rs654   A       G       21.4    PASS    DP=26   GT:X:Y:Z  1/0:0:1:0
chr3    10001   .       T       C       23.2    PASS    DP=29   GT:X:Y:Z  0/1:0:1:0

关键规则 #

• 若某变异与子集里两个变异的距离相同，可随机选择其中一个的元数据。
• 子集最后一个变异之后的所有全量变异，均使用该最后变异的元数据。
• 多样本场景下，每个样本的元数据仅从自身在子集里的记录获取，不能跨样本复用。

性能要求 #

数据集规模庞大（数百万变异、数千样本），需要内存友好、效率最优的实现方案。

已尝试的思路 #

优先考虑用bcftools实现，但不确定具体操作。曾尝试用以下命令提取元数据：

bcftools query -f '%CHROM %POS[\t%GT\t%X\t%Y\t%Z]\n' file.vcf > file_meta.txt

计划后续用Python处理，但未写出对应逻辑。

内容的提问来源于stack exchange，提问作者binf-er