合并 2 个匹配多个列的大文件并保留顺序(打印匹配和不匹配的值) - 从 awk 扩展

合并 2 个匹配多个列的大文件并保留顺序(打印匹配和不匹配的值) - 从 awk 扩展

我在合并 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 的快速内存,这也需要几个小时。

所以问题的关键是如何存储查找表。我建议使用二叉树来成倍减少执行时间。您可以使用perlC或其他某种语言来执行此操作,但此时此刻它会变得偏离主题。

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]
        }   
    }
}

相关内容