我在合并 2 个文件中的数据时遇到问题。它是带有染色体、位置、参考和替代等位基因的遗传数据,我需要通过匹配所有 4 列来合并文件 - 与参考和替代等位基因无论哪种方式。因此,我需要查找文件中的列 $1、$2、$4 和 $5 或 $1、$2、$5 和 $4 与数据文件中的列 $1、$5、$6 和 $7 完全匹配。在数据文件中保持确切的顺序非常重要 - 我无法对其进行排序(很遗憾不能使用 join - 这是我在此类问题的其他实例中找到的建议答案)。
我已经使用 awk 并让代码适用于具有几千行的示例文件,但它不能很好地扩展我的大型数据集(查找文件有 > 3 亿行;数据文件有 3000 万行) - 大概是因为代码需要保留内存 2 个用于查找的巨大数组。非常感谢您对可扩展代码(?在 perl 中?)的任何建议!
查找文件的格式为:
1 10150 rs371194064 C T
1 10165 rs796884232 A AC
1 10177 rs367896724 A AC
1 10177 rs201752861 A C
1 10180 rs201694901 T C
1 10199 rs905327004 A T
1 10231 rs200279319 C A
1 10234 rs145599635 C T
1 10235 rs540431307 T TA
1 10235 rs1035249121 T A
1 10235 rs1035249121 T C
1 10241 rs960927773 T C
1 10247 rs796996180 T C
1 10248 rs148908337 A T
1 10249 rs774211241 AAC A
我的数据文件的格式是
1 chr1 chr1:10177 1:10177_A_AC 10177 A AC
1 chr1 chr1:10235 1:10235_T_TA 10235 T TA
1 chr1 chr1:10352 1:10352_T_TA 10352 T TA
1 chr1 chr1:10505 1:10505_A_T 10505 A T
1 chr1 chr1:10506 1:10506_C_G 10506 C G
1 chr1 chr1:10511 1:10511_G_A 10511 G A
1 chr1 chr1:10539 1:10539_C_A 10539 C A
1 chr1 chr1:10542 1:10542_C_T 10542 C T
1 chr1 chr1:10579 1:10579_C_A 10579 C A
输出应如下所示:
1 rs367896724 1:10177_A_AC 10177 A AC A AC
1 rs540431307 1:10235_T_TA 10235 T TA T TA
1 chr1:10352 1:10352_T_TA 10352 T TA T TA
1 chr1:10505 1:10505_A_T 10505 A T A T
1 chr1:10506 1:10506_C_G 10506 C G C G
1 chr1:10511 1:10511_G_A 10511 G A G A
1 chr1:10539 1:10539_C_A 10539 C A C A
1 chr1:10542 1:10542_C_T 10542 C T C T
1 chr1:10579 1:10579_C_A 10579 C A C A
我设法为示例文件工作的 awk 代码如下:
awk 'BEGIN {OFS = "\t"} NR==FNR { #lookup file (323 million rows) key = $1 "," $2 "," $4 "," $5 present[key] = 1 ID[key] = $3 Ref[key] = $4 Alt[key] = $5 key1 = $1 "," $2 "," $4 "," $5 present1[key1] = 1 ID1[key1] = $3 Ref1[key1] = $4 Alt1[key1] = $5 next } { # my data file (3 million rows) key = $1 "," $5 "," $6 "," $7 key1 = $1 "," $5 "," $7 "," $6 if (present[key]) print $1, ID[key], $4, $5, $6, $7, Ref[key], Alt[key]; else if (present1[key1]) print $1, ID1[key1], $4, $5, $6, $7, Ref1[key1], Alt1[key1]; else print $1, $3, $4, $5, $6, $7, $6, $7 }' $lookupfile $mydatafile > $outputfile
答案1
问题不在于将该文件保留在内存中,而是扫描查找表以查找数据文件的每一行。您的代码没有显示它,但在幕后您执行了 3'000'000 次 323'000'000/2 = 几乎半千万亿字符串比较,在内存总线上移动了数千 TB。即使对于 200 GBit/s 的快速内存,这也需要几个小时。
所以问题的关键是如何存储查找表。我建议使用二叉树来成倍减少执行时间。您可以使用perl
或C
或其他某种语言来执行此操作,但此时此刻它会变得偏离主题。
unix 命令工具集无法帮助您解决此问题。
答案2
如果我的假设是正确的,那么两个文件都按染色体、碱基对位置、rs 编号(仅查找表)和最后的等位基因排序 - 至少显示的部分遵循此模式。在这种情况下,您不需要将整个查找表保留在内存中。相反,您只需要浏览每个文件一次,内存需求可以忽略不计:
依次遍历数据文件中的每个标记,然后在查找文件中搜索,直到找到匹配项或超出候选位置并确定没有匹配项。如果找到匹配项,则从查找表中提取相应的 rs 编号,否则仅使用数据表中的当前 chr:bp 组合。
使用下面的脚本我得到了您想要的输出。保存脚本,然后像这样使用它:
gawk -f scriptname datafile lookuptable outputfile
一些小的补充:为了获得有关处理的数据量的最小反馈,使用“#”和“.”。分别是数据表和查找表中每 10,000 行的输出。
#!/usr/bin/gawk -f
BEGIN {
OFS = "\t"
step = 10000
while (1==1) {
if ((getline indata < ARGV[1]) < 1)
break
if (!(na++ % step))
printf "\n#"
split(indata,a)
allequal = 0
while (1==1) {
if (!overrun) {
if ((getline inlookup < ARGV[2]) < 1)
break
if (!(nb++ % step))
printf "."
} else {
overrun=0
}
split(inlookup,b)
if (b[1]>a[1] || b[2]>a[5]) {
overrun=1
break
}
if (a[1]==b[1] && a[5]==b[2] && ((a[6]==b[4] && a[7]==b[5]) || (a[7]==b[4] && a[6]==b[5]))) {
allequal=1
break
}
}
if (allequal) {
print a[1],b[3],a[4],a[5],a[6],a[7],b[4],b[5] > ARGV[3]
} else {
print a[1],a[3],a[4],a[5],a[6],a[7],a[6],a[7] > ARGV[3]
}
}
}