linux 在另一个文件中查找大文件中的行的最快方法

hfsqlsce  于 2023-01-04  发布在  Linux
关注(0)|答案(6)|浏览(221)

我在while循环中使用grep来查找一个文件到另一个文件中的行,并将输出保存到一个新文件中。我的文件相当大(2.26亿行),脚本花费了很长时间(12天,还在继续)。你有什么建议来加快它吗?也许有比grep更好的方法?
(我还需要前面的行作为输出,因此grep-B 1。
下面是我的代码:

#!/bin/bash

while IFS= read -r line; do
  grep -B 1 $line K33.21mercounts.bf.trimmedreads.dumps.fa >> 21mercounts.bf.trimmedreads.diff.kmers.K33;
done <21mercounts.bf.trimmedreads.diff.kmers

更新:
包含要查找的行的输入文件大小为4.7 GB,行数为2.26亿,如下所示:

AAAGAAAAAAAAAGCTAAAAT
ATCTCGACGCTCATCTCAGCA
GTTCGTCGGAGAGGAGAGAAC
GAGGACTATAAAATTGTCGCA
GGCTTCAATAATTTGTATAAC
GACATAGAATCACGAGTGACC
TGGTGAGTGACATCCTTGACA
ATGAAAACTGCCAGCAAACTC
AAAAAACTTACCTTAAAAAGT
TTAGTACACAATATCTCCCAA

要查找的文件大小为26 GB,包含20亿行,如下所示:

>264638
AAAAAAAAAAAAAAAAAAAAA
>1
AAAGAAAAAAAAAGCTAAAAT
>1
ATCTCGACGCTCATCTCAGCA
>1
GTTCGTCGGAGAGGAGAGAAC
>28
TCTTTTCAGGAGTAATAACAA
>13
AATCATTTTCCGCTGGAGAGA
>38
ATTCAATAAATAATAAATTAA
>2
GAGGACTATAAAATTGTCGCA
>1
GGCTTCAATAATTTGTATAAC

预期输出如下:

>1
AAAGAAAAAAAAAGCTAAAAT
>1
ATCTCGACGCTCATCTCAGCA
>1
GTTCGTCGGAGAGGAGAGAAC
>2
GAGGACTATAAAATTGTCGCA
>1
GGCTTCAATAATTTGTATAAC
waxmsbnn

waxmsbnn1#

您可以尝试不使用shell循环并使用固定字符串搜索的grep -f命令:

grep -B1 -Ff 21mercounts.bf.trimmedreads.diff.kmers \
 K33.21mercounts.bf.trimmedreads.dumps.fa > 21mercounts.bf.trimmedreads.diff.kmers.K33
jc3wubiy

jc3wubiy2#

如果不需要保留原始顺序,使用GNU uniq和GNU sed

{ cat 21mercounts.bf.trimmedreads.diff.kmers
  sed -n 'x;n;G;s/\n//p' K33.21mercounts.bf.trimmedreads.dumps.fa
} | LC_ALL=C sort | uniq -w21 -D |
sed -n 's/\(.*\)>\(.*\)/>\2\n\1/p' > 21mercounts.bf.trimmedreads.diff.kmers.K33
llycmphe

llycmphe3#

有不少工具(例如ripgrep)和选项(-f-F-x)来加速您的基本方法。但它们基本上都是与您现在使用的方法一样慢的方法,"只是"加速了一个巨大但仍然 * 恒定 * 的因子。对于您的问题和输入大小,我建议完全改变这种方法。有很多不同的方法可以解决你的问题。首先,让我们定义一些变量来估计这些方法的加速比:

问题

一个26 GBhaystack文件,其中h= 1百万个条目(描述、序列)= 20亿行,例如

>28
TCTTTTCAGGAGTAATAACAA
>13
AATCATTTTCCGCTGGAGAGA
>38
ATTCAATAAATAATAAATTAA
...

4.7GB文件,n= 2.26亿行,每行长度m= 21,例如

GACATAGAATCACGAGTGACC
TGGTGAGTGACATCCTTGACA
ATGAAAACTGCCAGCAAACTC
...

对于所有的needle,我们希望从干草堆中提取相应的条目(如果它们存在的话)。

溶液

我们假设n〈h且m为常数,则O(n + h)= O(h),O(m)= O(1),依此类推。
我们的目标是最小化迭代最大文件(=干草堆)的次数。

未处理-O(h·n)时间

目前,您使用的是简单的方法,对于每根针,整个干草堆都搜索一次。

将针放入数据结构;搜索一次干草堆-O(h)时间

将所有指针存储在一个具有快速contains()操作的数据结构中,然后迭代干草堆并为每个条目调用needles.contains(haystackEntry),以确定它是否是您要搜索的内容。
目前,您的"数据结构"是一个列表,"构建"需要O(1)时间(因为它已经是这种形式),但是查询 * 一次 * 需要O(n)时间!
下面的数据结构需要O(n)的时间来填充,O(1)的时间来查询一次,导致O(n + h·1)= O(h)的时间。

  • Trie s(=前缀树)可以表示为正则表达式,所以你可以坚持使用grep,例如指针ABCABXXBC可以存储在Trie正则表达式^(AB(C|X)|XBC)中,但是在bash中将指针列表转换为这样的Trie正则表达式有点复杂。
  • 散列Map在awk中可用,请参见sundeep's answer。但是将4.7 GB的原始数据放在这样的内存结构中可能效率不高(取决于可用内存。散列Map需要比原始数据大很多倍)。

无论哪种方式,数据结构和bash都不能很好地混合。即使我们切换到一种更好的语言,我们也必须在每次程序运行时重新构建或存储/加载结构。

排序一切;搜索一次干草堆-O(h·log(h)+h)时间

首先对干草堆和针进行排序,然后只迭代干草堆一次。
用第一根针从头开始搜索干草堆。当到达一个干草堆条目时,必须在当前针后面排序,用下一根针从当前位置继续搜索。
这在bash中很容易做到,这里我们使用GNU coreutils来使处理更容易、更快、更安全:

export LC_ALL=C  # speeds up sorting
mem=66%    # Max. memory to be used while sorting. More is better.
sep=$'\f'  # A character not appearing in your data.

paste -d"$sep" - -  < haystack > haystack2

sort -S66% -o needles2 needles
sort -t"$sep" -k2,2 -S"$mem" -o haystack2 haystack2

# --nocheck-order is not needed, but speeds up the process
join -t"$sep" -22 -o2.1,2.2 --nocheck-order needles2 haystack2 |
tr "$sep" \\n

这将改变输出的顺序。如果需要原始顺序的输出,请使用Schwartzian transform(= decorate-sort-undecorate):在对针/干草堆进行排序之前,存储它们的行号。在整个过程中拖动这些行号。最后,按行号对找到的条目进行排序。最后,删除行号并打印结果。

export LC_ALL=C  # speeds up sorting
mem=66%    # Max. memory to be used while sorting. More is better.
sep=$'\f'  # A character not appearing in your data.

nl -ba -d '' -s"$sep" needles > needles2
paste -d"$sep" - -  < haystack | nl -ba -d '' -s"$sep" > haystack2

sort -t"$sep" -k2,2 -S"$mem" -o needles2 needles2
sort -t"$sep" -k3,3 -S"$mem" -o haystack2 haystack2

# --nocheck-order is not needed, but speeds up the process
join -t"$sep" -12 -23 -o1.1,2.1,2.2,2.3 --nocheck-order needles2 haystack2 > result
sort -t"$sep" -k1,2n -S"$mem" -o result result
cut -d"$sep" -f3- result | tr "$sep" \\n
js81xvg6

js81xvg64#

这里有一个使用awk的解决方案。不确定它是否会比grepripgrep更快,但由于基于哈希的查找,这是可能的。这里假设您的RAM足够大,可以加载第一个文件(4.7 GB和2.26亿行)。

$ awk 'NR==FNR{a[$1]; next} $0 in a{print p; print} {p=$0}' f1 f2
>1
AAAGAAAAAAAAAGCTAAAAT
>1
ATCTCGACGCTCATCTCAGCA
>1
GTTCGTCGGAGAGGAGAGAAC
>2
GAGGACTATAAAATTGTCGCA
>1
GGCTTCAATAATTTGTATAAC

mawk通常是最快的选项,但我遇到过gawk更快的例子,特别是对于像下面这样的阵列。如果你可以安装frawk,那会给予你更快的结果。命令需要稍微修改一下:

frawk 'NR==FNR{a[$1]; next} $0 in a{print p; print $0} {p=$0}' f1 f2
t98cgbkg

t98cgbkg5#

每当我处理这么大的文件时,我几乎总是要对它们进行排序。排序很慢,但比while read循环扫描20亿行2.26亿次要少得多。

sort 4GB>4gb.srt

以及

sed '/>/{N;s/\n/ /}' 26GB |sort -t' ' -k2 >25gb.srt

它将生成一个如下所示的文件:

>264638 AAAAAAAAAAAAAAAAAAAAA
>1 AAAGAAAAAAAAAGCTAAAAT
>13 AATCATTTTCCGCTGGAGAGA
>1 ATCTCGACGCTCATCTCAGCA
>38 ATTCAATAAATAATAAATTAA
>2 GAGGACTATAAAATTGTCGCA
>1 GGCTTCAATAATTTGTATAAC
>1 GTTCGTCGGAGAGGAGAGAAC
>28 TCTTTTCAGGAGTAATAACAA

现在您只需通读每个文件一次。

$ cat tst
awk 'BEGIN{ getline key < "4gb.srt"; }
 $2  < key { next; }
 $2  > key { while ($2 > key){ getline key < "4gb.srt"; } }
 $2 == key {  $0=gensub(/ /,"\n",1); print }' 25gb.srt

$ ./tst
>1
AAAGAAAAAAAAAGCTAAAAT
>1
ATCTCGACGCTCATCTCAGCA
>2
GAGGACTATAAAATTGTCGCA
>1
GGCTTCAATAATTTGTATAAC
>1
GTTCGTCGGAGAGGAGAGAAC

排序和你的不一样,但除此之外那还管用吗?
(Try先使用较小的文件进行一些测试...)

附录

请参考Socowi's better implementation,但我被要求解释awk,所以-
首先,请看上面,我将较大的“haystraw”文件解析为按关键字字段排序的单行,这将是$2 in awk,并将较小的“needles”文件解析为相同的顺序。

BEGIN{ getline key < "4gb.srt"; }

这只是通过阅读适当的文件将第一个“needle”初始化为一个名为key的变量。
然后,当awk读取“haystraw”文件的每一行时,它自动将其解析为字段--因为我们将它们堆叠起来,所以第一个字段是原始干草堆的前一行,第二个字段是要检查的值,因此我们在key$2之间进行比较。

$2  < key { next; } # skip ahead to next key/needle

如果当前吸管小于针头,则将其扔掉并抓住下一根。

$2  > key { while ($2 > key){ getline key < "4gb.srt"; } }

如果当前吸管比针“大”,那么针不在文件中,下一个吸管也可能不在文件中,所以我们按顺序抓取针,然后进行比较,直到它们赶上为止。
这里实际上有一个潜在的bug--它无法确认某些内容已被读取,并且可能在指针用完时挂在一个无限循环中。这一节 * 应该 * 是这样的--

$2  > key { while ( ($2 > key) { if( 0 == getline key < "4gb.srt" ) key = "ZZZZZZZZZZZZZZZZZZZZZZ"; } }

最后,

$2 == key {  $0=gensub(/ /,"\n",1); print }' 25gb.srt

如果它们匹配,则在前一记录和匹配行之间重新插入换行符,并打印它们。
真的也应该有一个END{ close("4gb.srt") }

ddrv8njm

ddrv8njm6#

grep可以同时搜索多个模式(在单独的文件中给出),因此阅读K33.21mercounts.bf.trimmedreads.dumps.fa只需要执行一次。

#!/bin/bash

grep --f 21mercounts.bf.trimmedreads.diff.kmers -B 1 K33.21mercounts.bf.trimmedreads.dumps.fa >> 21mercounts.bf.trimmedreads.diff.kmers.K33;

但是,它可能需要大量RAM

相关问题