计算数百 GB 数据中的子序列
Count subsequences in hundreds of GB of data
我正在尝试处理一个非常大的文件并计算文件中特定长度的所有序列的频率。
为了说明我在做什么，考虑一个包含序列 abcdefabcgbacbdebdbbcaebfebfebfeb
的小输入文件
下面，代码读入整个文件，取长度为n的第一个子串（下面我把这个设置为5，虽然我想能够改变这个）并计算它的频率：
abcde => 1
下一行，它向右移动一个字符并做同样的事情：
bcdef => 1
然后继续处理字符串的其余部分并打印 5 个最频繁的序列：
open my $in, '<', 'in.txt' or die $!; # 'abcdefabcgbacbdebdbbcaebfebfebfeb'
my $seq = <$in>; # read whole file into string
my $len = length($seq);
my $seq_length = 5; # set k-mer length
my %data;
for (my $i = 0; $i <= $len - $seq_length; $i++) {
my $kmer = substr($seq, $i, $seq_length);
$data{$kmer}++;
}
# print the hash, showing only the 5 most frequent k-mers
my $count = 0;
foreach my $kmer (sort { $data{$b} <=> $data{$a} } keys %data ){
print "$kmer $data{$kmer}\n";
$count++;
last if $count >= 5;
}
ebfeb 3
febfe 2
bfebf 2
bcaeb 1
abcgb 1
不过，我想找到一种更有效的方法来实现这一点。如果输入文件是 10GB 或 1000GB，那么将整个文件读入一个字符串将非常耗费内存。
我考虑过读取字符块，比如一次读取 100 个字符，然后按上述方法继续，但在这里，跨越 2 个块的序列将无法正确计数。
我的想法是只从字符串中读取 n 个字符，然后移动到下一个 n 个字符并执行相同的操作，如上所示在哈希中计算它们的频率。
- 关于我如何做到这一点有什么建议吗？我看过 read 使用偏移量，但无法理解如何将其合并到此处
substr 是这项任务最有效的内存工具吗？
从您自己的代码来看，您的数据文件看起来只有一行数据——没有被换行符分隔——所以我在下面的解决方案中假设了这一点。即使该行末尾可能有一个换行符，最后选择五个最频繁的子序列也会将其排除在外，因为它只发生一次
该程序使用sysread从文件中获取任意大小的数据块并将其附加到我们内存中已有的数据中
循环体大部分与您自己的代码相似，但我使用了 for 的列表版本而不是 C 风格的，因为它更清晰
处理完每个块后，内存中的数据被截断为最后 SEQ_LENGTH-1 个字节，然后循环的下一个周期从文件中提取更多数据
我还对 K-mer 大小和块大小使用了常量。毕竟它们是不变的！
输出数据是在CHUNK_SIZE设置为7的情况下产生的，这样会有很多跨界子序列的实例。它匹配您自己所需的输出，但最后两个条目的计数为 1。这是因为 Perl 的散列键固有的随机顺序，如果您需要具有相同计数的特定序列顺序，则必须指定它，以便我可以更改排序
use strict;
use warnings 'all';
use constant SEQ_LENGTH => 5; # K-mer length
use constant CHUNK_SIZE => 1024 * 1024; # Chunk size - say 1MB
my $in_file = shift // 'in.txt';
open my $in_fh, '<', $in_file or die qq{Unable to open "$in_file" for input: $!};
my %data;
my $chunk;
my $length = 0;
while ( my $size = sysread $in_fh, $chunk, CHUNK_SIZE, $length ) {
$length += $size;
for my $offset ( 0 .. $length - SEQ_LENGTH ) {
my $kmer = substr $chunk, $offset, SEQ_LENGTH;
++$data{$kmer};
}
$chunk = substr $chunk, -(SEQ_LENGTH-1);
$length = length $chunk;
}
my @kmers = sort { $data{$b} <=> $data{$a} } keys %data;
print "$_ $data{$_}\n" for @kmers[0..4];
输出
ebfeb 3
febfe 2
bfebf 2
gbacb 1
acbde 1
注意行：$chunk = substr $chunk, -(SEQ_LENGTH-1); 当我们通过 while 循环时设置 $chunk。这可确保正确计算跨越 2 个块的字符串。
$chunk = substr $chunk, -4 语句从当前块中删除除最后四个字符以外的所有字符，以便下一次读取将文件中的 CHUNK_SIZE 字节附加到剩余的字符。这样搜索将继续，但除了下一个块之外，还会从前一个块的最后 4 个字符开始：数据不会落入块之间的 "crack"。
最直接的方法是使用substr()函数：
% time perl -e '$/ = 48576;
while ($s = <>) { for $i (0..length $s) {
$hash{ substr($s, $i, 5) }++ } }
foreach my $k (sort { $hash{$b} <=> $hash{$a} } keys %hash) {
print "$k $hash{$k}\n"; $it++; last if $it == 5;}' nucleotide.data
NNCTA 337530
GNGGA 337362
NCACT 337304
GANGN 337290
ACGGC 337210
269.79 real 268.92 user 0.66 sys
iterating along a string 上的 Perl Monks 节点是有用的资源，@Jonathan Leffler、@ÆvarArnfjörðBjarmason、@Vorsprung、@ThisSuitIsBlackNotm @ 的回复和评论也是有用的资源borodin 和@ikegami 在这个 SO posting 中。正如所指出的，非常大的文件的问题是内存，这反过来又需要以块的形式读取文件。当以块的形式从文件中读取时，如果您的代码正在遍历数据，它必须正确处理从一个 chunk/source 到下一个的切换，而不会丢失任何字节。
作为一个简单的例子，next unless length $kmer == 5; 将在上面脚本中的每个 1048576 byte/character 迭代期间被检查，这意味着字符串存在于一个块的末尾和块的开头另一个将被遗漏（cf. @ikegami 和@Borodin 的解决方案）。这将改变结果计数，尽管可能不会以统计显着的方式[1]。 @borodin 和@ikegami 都解决了块之间 missing/overlapping 字符串的问题，方法是将每个块附加到前一个块的剩余字符，因为它们在 while() 循环中 sysread 。请参阅 Borodin 的回复和评论以了解其工作原理。
使用Stream::Reader
既然 perl 已经存在了很长一段时间并且收集了很多有用的代码，另一个完全有效的方法是寻找一个实现相同目的的 CPAN 模块。 Stream::Reader 可以创建一个 "stream" 文件句柄接口，将分块问题的解决方案包装在一组方便的数据访问函数后面。
use Stream::Reader;
use strict;
use warnings;
open( my $handler, "<", shift );
my $stream = Stream::Reader->new( $handler, { Mode => "UB" } );
my %hash;
my $string;
while ($stream->readto("\n", { Out => $string }) ) {
foreach my $i (0..length $string) {
$hash{ substr($string, $i, 5) }++
}
}
my $it;
foreach my $k (sort { $hash{$b} <=> $hash{$a} } keys %hash ) {
print "$k $hash{$k}\n";
$it++; last if $it == 5;
}
在测试数据文件 nucleotide.data 上，Borodin 的脚本和上面显示的 Stream::Reader 方法都产生了相同的前五名结果。请注意与上面 shell 命令的结果相比的细微差别。这说明了正确处理块读取数据的必要性。
NNCTA 337530
GNGGA 337362
NCACT 337305
GANGN 337290
ACGGC 337210
基于 Stream::Reader 的脚本明显更快：
time perl sequence_search_stream-reader.pl nucleotide.data
252.12s
time perl sequence_search_borodin.pl nucleotide.data
350.57s
文件 nucleotide.data 大小为 1Gb，由大约 10 亿个字符的单个字符串组成：
% wc nucleotide.data
0 0 1048576000 nucleotide.data
% echo `head -c 20 nucleotide.data`
NCCANGCTNGGNCGNNANNA
我使用这个命令来创建文件：
perl -MString::Random=random_regex -e '
open (my $fh, ">>", "nucleotide.data");
for (0..999) { print $fh random_regex(q|[GCNTA]{1048576}|) ;}'
列表和字符串
由于应用程序应该一次读取一个块并沿着数据的长度移动这个 $seq_length 大小的 window 构建一个用于跟踪字符串频率的散列，我认为 "lazy list" 方法可能在这里起作用。但是，要通过一组数据移动 window（或 slide as with List::Gen) reading elements natatime，需要一个列表。
我将数据视为一个非常长的字符串，首先必须将其制成列表才能使这种方法起作用。我不确定这样做的效率如何。不过，这是我尝试 "lazy list" 解决问题的方法：
use List::Gen 'slide';
$/ = 48575; # Read a million character/bytes at a time.
my %hash;
while (my $seq = <>) {
chomp $seq;
foreach my $kmer (slide { join("", @_) } 5 => split //, $seq) {
next unless length $kmer == 5;
$hash{$kmer}++;
}
}
foreach my $k (sort { $hash{$b} <=> $hash{$a} } keys %hash) {
print "$k $hash{$k}\n";
$it++; last if $it == 5;
}
我不确定这是 "typical perl"（当然是 TIMTOWDI），我想还有其他技术（cf. gather/take）和实用程序适合这个任务。我最喜欢@Borodin 的回复，因为它似乎是执行此任务的最常见方式，并且对于提到的可能较大的文件大小 (100Gb) 而言效率更高。
是否有 fast/best 方法将字符串转换为列表或对象？使用增量 read() 或 sysread() 和 substr 在这一点上获胜，但即使使用 sysread 一个 1000Gb 的字符串也需要大量内存来存储结果哈希。或许 serialized/cached 将散列存储到磁盘的技术随着它的增长超过一定大小将适用于非常非常大的字符串，这些字符串很容易创建非常大的散列。
后记和结果
List::Gen 方法始终比@Borodin 的方法慢 5 到 6 倍。最快的脚本使用了 Stream::Reader 模块。结果是一致的，每个脚本都选择了相同的前五个字符串和两个较小的文件：
100万字符核苷酸串
sequence_search_stream-reader.pl : 0.26s
sequence_search_borodin.pl : 0.39s
sequence_search_listgen.pl : 2.04s
8300万字符核苷酸串
与文件xaa中的数据：
wc xaa
0 1 83886080 xaa
% time perl sequence_search_stream-reader.pl xaa
GGCNG 31510
TAGNN 31182
AACTA 30944
GTCAN 30792
ANTAT 30756
21.33 real 20.95 user 0.35 sys
% time perl sequence_search_borodin.pl xaa
GGCNG 31510
TAGNN 31182
AACTA 30944
GTCAN 30792
ANTAT 30756
28.13 real 28.08 user 0.03 sys
% time perl sequence_search_listgen.pl xaa
GGCNG 31510
TAGNN 31182
AACTA 30944
GTCAN 30792
ANTAT 30756
157.54 real 156.93 user 0.45 sys
10亿字符核苷酸串
在更大的文件中，差异幅度相似，但是，因为它没有正确处理跨越块边界的序列，所以 List::Gen 脚本与 shell 命令行有相同的差异在此 post 的开头。较大的文件意味着许多块边界和计数差异。
sequence_search_stream-reader.pl : 252.12s
sequence_search_borodin.pl : 350.57s
sequence_search_listgen.pl : 1928.34s
块边界问题当然可以解决，但我很想知道使用 "lazy list" 方法引入的其他潜在错误或瓶颈。如果使用 slide 到 "lazily" 沿字符串移动在 CPU 用法方面有任何好处，它似乎因为需要在字符串中列出列表而变得毫无意义开始。
我对跨块边界读取数据留作实施练习并不感到惊讶（也许它无法处理 "magically"），但我想知道可能存在哪些其他 CPAN 模块或陈旧的子例程样式解决方案。
1. 在读取 TB 文件的每兆字节末尾跳过四个字符 - 因此有四个 5 字符的字符串组合意味着结果将不包括 3/10000 来自最终计数的 1%。
echo "scale=10; 100 * (1024^4/1024^2 ) * 4 / 1024^4 " | bc
.0003814697
一般来说，Perl 在逐字符处理解决方案（如上面发布的那些解决方案）方面确实很慢，但在诸如正则表达式之类的解决方案方面快得多，因为基本上你的开销主要是多少您正在执行的运算符。
因此，如果您可以将其转变为基于正则表达式的解决方案，那就更好了。
尝试这样做：
$ perl -wE 'my $str = "abcdefabcgbacbdebdbbcaebfebfebfeb"; for my $pos (0..4) { $str =~ s/^.// if $pos; say for $str =~ m/(.{5})/g }'|sort|uniq -c|sort -nr|head -n 5
3 ebfeb
2 febfe
2 bfebf
1 gbacb
1 fabcg
即我们在 $str 中有我们的字符串，然后我们传递它 5 次生成 5 个字符的序列，在第一次传递之后我们开始从字符串的前面切掉一个字符。在很多语言中，这会 真的 慢，因为你必须重新分配整个字符串，但是 perl 会针对这种特殊情况作弊，只是将字符串的索引设置为 1 + 之前是什么。
我没有对此进行基准测试，但我敢打赌，这样做比上面的算法更可行，你也可以在 perl 中进行 uniq 计数，当然是通过递增哈希（使用 / e regex 选项可能是最快的方式），但我只是将其卸载到此实现中的 |sort|uniq -c，这可能更快。
在 perl 中完成这一切的稍微改变的实现：
$ perl -wE 'my $str = "abcdefabcgbacbdebdbbcaebfebfebfeb"; my %occur; for my $pos (0..4) { substr($str, 0, 1) = "" if $pos; $occur{$_}++ for $str =~ m/(.{5})/gs }; for my $k (sort { $occur{$b} <=> $occur{$a} } keys %occur) { say "$occur{$k} $k" }'
3 ebfeb
2 bfebf
2 febfe
1 caebf
1 cgbac
1 bdbbc
1 acbde
1 efabc
1 aebfe
1 ebdbb
1 fabcg
1 bacbd
1 bcdef
1 cbdeb
1 defab
1 debdb
1 gbacb
1 bdebd
1 cdefa
1 bbcae
1 bcgba
1 bcaeb
1 abcgb
1 abcde
1 dbbca
其背后代码的漂亮格式：
my $str = "abcdefabcgbacbdebdbbcaebfebfebfeb";
my %occur;
for my $pos (0..4) {
substr($str, 0, 1) = "" if $pos;
$occur{$_}++ for $str =~ m/(.{5})/gs;
}
for my $k (sort { $occur{$b} <=> $occur{$a} } keys %occur) {
say "$occur{$k} $k";
}
即使您在处理文件之前没有将整个文件读入内存，您仍然可能 运行 内存不足。
一个 10 GiB 的文件包含将近 11E9 个序列。
如果您的序列是从一组 5 个字符中选出的 5 个字符的序列，则只有 55 = 3,125 个唯一序列，这很容易存储在内存中。
如果您的序列是从一组 5 个字符中选出的 20 个字符的序列，则有 520 = 95E12 个唯一序列，因此一个 10 GiB 文件的所有 11E9 个序列可以独一无二。那不适合记忆。
在这种情况下，我建议执行以下操作：
创建一个包含原始文件所有序列的文件。
以下是分块读取文件，而不是一次全部读取。棘手的部分是处理跨越两个块的序列。以下程序使用 sysread[1] 从文件中获取任意大小的数据块并将其附加到先前读取块的最后几个字符。最后一个细节允许对跨越块的序列进行计数。
perl -e'
use strict;
use warnings qw( all );
use constant SEQ_LENGTH => 20;
use constant CHUNK_SIZE => 1024 * 1024;
my $buf = "";
while (1) {
my $size = sysread(\*STDIN, $buf, CHUNK_SIZE, length($buf));
die($!) if !defined($size);
last if !$size;
for my $offset ( 0 .. length($buf) - SEQ_LENGTH ) {
print(substr($buf, $offset, SEQ_LENGTH), "\n");
}
substr($buf, 0, -(SEQ_LENGTH-1), "");
}
' <in.txt >sequences.txt
对序列进行排序。
sort sequences.txt >sorted_sequences.txt
计算每个序列的实例数，并将计数与序列一起存储在另一个文件中。
perl -e'
use strict;
use warnings qw( all );
my $last = "";
my $count;
while (<>) {
chomp;
if ($_ eq $last) {
++$count;
} else {
print("$count $last\n") if $count;
$last = $_;
$count = 1;
}
}
' sorted_sequences.txt >counted_sequences.txt
按计数对序列排序。
sort -rns counted_sequences.txt >sorted_counted_sequences.txt
提取结果。
perl -e'
use strict;
use warnings qw( all );
my $last_count;
while (<>) {
my ($count, $seq) = split;
last if $. > 5 && $count != $last_count;
print("$seq $count\n");
$last_count = $count;
}
' sorted_counted_sequences.txt
这也打印出并列第 5 名。
这可以通过调整传递给 sort[2] 的参数来优化，但它应该提供不错的性能。
sysread 比之前建议的 read 更快，因为后者在内部执行一系列 4 KiB 或 8 KiB 读取（取决于您的 Perl 版本）。
鉴于序列的固定长度性质，您还可以将序列压缩为 ceil(log256(520)) = 6 个字节，然后将它们进行 base64 编码为 ceil(6 * 4/3) = 8 个字节。这意味着每个序列需要的字节减少 12 个，大大减少了读写量。
此答案的部分内容改编自 user:622310 根据 cc by-sa 3.0 许可的内容。
我正在尝试处理一个非常大的文件并计算文件中特定长度的所有序列的频率。
为了说明我在做什么，考虑一个包含序列 abcdefabcgbacbdebdbbcaebfebfebfeb
下面，代码读入整个文件，取长度为n的第一个子串（下面我把这个设置为5，虽然我想能够改变这个）并计算它的频率：
abcde => 1
下一行，它向右移动一个字符并做同样的事情：
bcdef => 1
然后继续处理字符串的其余部分并打印 5 个最频繁的序列：
open my $in, '<', 'in.txt' or die $!; # 'abcdefabcgbacbdebdbbcaebfebfebfeb'
my $seq = <$in>; # read whole file into string
my $len = length($seq);
my $seq_length = 5; # set k-mer length
my %data;
for (my $i = 0; $i <= $len - $seq_length; $i++) {
my $kmer = substr($seq, $i, $seq_length);
$data{$kmer}++;
}
# print the hash, showing only the 5 most frequent k-mers
my $count = 0;
foreach my $kmer (sort { $data{$b} <=> $data{$a} } keys %data ){
print "$kmer $data{$kmer}\n";
$count++;
last if $count >= 5;
}
ebfeb 3
febfe 2
bfebf 2
bcaeb 1
abcgb 1
不过，我想找到一种更有效的方法来实现这一点。如果输入文件是 10GB 或 1000GB，那么将整个文件读入一个字符串将非常耗费内存。
我考虑过读取字符块，比如一次读取 100 个字符，然后按上述方法继续，但在这里，跨越 2 个块的序列将无法正确计数。
我的想法是只从字符串中读取 n 个字符，然后移动到下一个 n 个字符并执行相同的操作，如上所示在哈希中计算它们的频率。
- 关于我如何做到这一点有什么建议吗？我看过 read 使用偏移量，但无法理解如何将其合并到此处
substr是这项任务最有效的内存工具吗？
从您自己的代码来看，您的数据文件看起来只有一行数据——没有被换行符分隔——所以我在下面的解决方案中假设了这一点。即使该行末尾可能有一个换行符，最后选择五个最频繁的子序列也会将其排除在外，因为它只发生一次
该程序使用sysread从文件中获取任意大小的数据块并将其附加到我们内存中已有的数据中
循环体大部分与您自己的代码相似，但我使用了 for 的列表版本而不是 C 风格的，因为它更清晰
处理完每个块后，内存中的数据被截断为最后 SEQ_LENGTH-1 个字节，然后循环的下一个周期从文件中提取更多数据
我还对 K-mer 大小和块大小使用了常量。毕竟它们是不变的！
输出数据是在CHUNK_SIZE设置为7的情况下产生的，这样会有很多跨界子序列的实例。它匹配您自己所需的输出，但最后两个条目的计数为 1。这是因为 Perl 的散列键固有的随机顺序，如果您需要具有相同计数的特定序列顺序，则必须指定它，以便我可以更改排序
use strict;
use warnings 'all';
use constant SEQ_LENGTH => 5; # K-mer length
use constant CHUNK_SIZE => 1024 * 1024; # Chunk size - say 1MB
my $in_file = shift // 'in.txt';
open my $in_fh, '<', $in_file or die qq{Unable to open "$in_file" for input: $!};
my %data;
my $chunk;
my $length = 0;
while ( my $size = sysread $in_fh, $chunk, CHUNK_SIZE, $length ) {
$length += $size;
for my $offset ( 0 .. $length - SEQ_LENGTH ) {
my $kmer = substr $chunk, $offset, SEQ_LENGTH;
++$data{$kmer};
}
$chunk = substr $chunk, -(SEQ_LENGTH-1);
$length = length $chunk;
}
my @kmers = sort { $data{$b} <=> $data{$a} } keys %data;
print "$_ $data{$_}\n" for @kmers[0..4];
输出
ebfeb 3
febfe 2
bfebf 2
gbacb 1
acbde 1
注意行：$chunk = substr $chunk, -(SEQ_LENGTH-1); 当我们通过 while 循环时设置 $chunk。这可确保正确计算跨越 2 个块的字符串。
$chunk = substr $chunk, -4 语句从当前块中删除除最后四个字符以外的所有字符，以便下一次读取将文件中的 CHUNK_SIZE 字节附加到剩余的字符。这样搜索将继续，但除了下一个块之外，还会从前一个块的最后 4 个字符开始：数据不会落入块之间的 "crack"。
最直接的方法是使用substr()函数：
% time perl -e '$/ = 48576;
while ($s = <>) { for $i (0..length $s) {
$hash{ substr($s, $i, 5) }++ } }
foreach my $k (sort { $hash{$b} <=> $hash{$a} } keys %hash) {
print "$k $hash{$k}\n"; $it++; last if $it == 5;}' nucleotide.data
NNCTA 337530
GNGGA 337362
NCACT 337304
GANGN 337290
ACGGC 337210
269.79 real 268.92 user 0.66 sys
iterating along a string 上的 Perl Monks 节点是有用的资源，@Jonathan Leffler、@ÆvarArnfjörðBjarmason、@Vorsprung、@ThisSuitIsBlackNotm @ 的回复和评论也是有用的资源borodin 和@ikegami 在这个 SO posting 中。正如所指出的，非常大的文件的问题是内存，这反过来又需要以块的形式读取文件。当以块的形式从文件中读取时，如果您的代码正在遍历数据，它必须正确处理从一个 chunk/source 到下一个的切换，而不会丢失任何字节。
作为一个简单的例子，next unless length $kmer == 5; 将在上面脚本中的每个 1048576 byte/character 迭代期间被检查，这意味着字符串存在于一个块的末尾和块的开头另一个将被遗漏（cf. @ikegami 和@Borodin 的解决方案）。这将改变结果计数，尽管可能不会以统计显着的方式[1]。 @borodin 和@ikegami 都解决了块之间 missing/overlapping 字符串的问题，方法是将每个块附加到前一个块的剩余字符，因为它们在 while() 循环中 sysread 。请参阅 Borodin 的回复和评论以了解其工作原理。
使用Stream::Reader
既然 perl 已经存在了很长一段时间并且收集了很多有用的代码，另一个完全有效的方法是寻找一个实现相同目的的 CPAN 模块。 Stream::Reader 可以创建一个 "stream" 文件句柄接口，将分块问题的解决方案包装在一组方便的数据访问函数后面。
use Stream::Reader;
use strict;
use warnings;
open( my $handler, "<", shift );
my $stream = Stream::Reader->new( $handler, { Mode => "UB" } );
my %hash;
my $string;
while ($stream->readto("\n", { Out => $string }) ) {
foreach my $i (0..length $string) {
$hash{ substr($string, $i, 5) }++
}
}
my $it;
foreach my $k (sort { $hash{$b} <=> $hash{$a} } keys %hash ) {
print "$k $hash{$k}\n";
$it++; last if $it == 5;
}
在测试数据文件 nucleotide.data 上，Borodin 的脚本和上面显示的 Stream::Reader 方法都产生了相同的前五名结果。请注意与上面 shell 命令的结果相比的细微差别。这说明了正确处理块读取数据的必要性。
NNCTA 337530
GNGGA 337362
NCACT 337305
GANGN 337290
ACGGC 337210
基于 Stream::Reader 的脚本明显更快：
time perl sequence_search_stream-reader.pl nucleotide.data
252.12s
time perl sequence_search_borodin.pl nucleotide.data
350.57s
文件 nucleotide.data 大小为 1Gb，由大约 10 亿个字符的单个字符串组成：
% wc nucleotide.data
0 0 1048576000 nucleotide.data
% echo `head -c 20 nucleotide.data`
NCCANGCTNGGNCGNNANNA
我使用这个命令来创建文件：
perl -MString::Random=random_regex -e '
open (my $fh, ">>", "nucleotide.data");
for (0..999) { print $fh random_regex(q|[GCNTA]{1048576}|) ;}'
列表和字符串
由于应用程序应该一次读取一个块并沿着数据的长度移动这个 $seq_length 大小的 window 构建一个用于跟踪字符串频率的散列，我认为 "lazy list" 方法可能在这里起作用。但是，要通过一组数据移动 window（或 slide as with List::Gen) reading elements natatime，需要一个列表。
我将数据视为一个非常长的字符串，首先必须将其制成列表才能使这种方法起作用。我不确定这样做的效率如何。不过，这是我尝试 "lazy list" 解决问题的方法：
use List::Gen 'slide';
$/ = 48575; # Read a million character/bytes at a time.
my %hash;
while (my $seq = <>) {
chomp $seq;
foreach my $kmer (slide { join("", @_) } 5 => split //, $seq) {
next unless length $kmer == 5;
$hash{$kmer}++;
}
}
foreach my $k (sort { $hash{$b} <=> $hash{$a} } keys %hash) {
print "$k $hash{$k}\n";
$it++; last if $it == 5;
}
我不确定这是 "typical perl"（当然是 TIMTOWDI），我想还有其他技术（cf. gather/take）和实用程序适合这个任务。我最喜欢@Borodin 的回复，因为它似乎是执行此任务的最常见方式，并且对于提到的可能较大的文件大小 (100Gb) 而言效率更高。
是否有 fast/best 方法将字符串转换为列表或对象？使用增量 read() 或 sysread() 和 substr 在这一点上获胜，但即使使用 sysread 一个 1000Gb 的字符串也需要大量内存来存储结果哈希。或许 serialized/cached 将散列存储到磁盘的技术随着它的增长超过一定大小将适用于非常非常大的字符串，这些字符串很容易创建非常大的散列。
后记和结果
List::Gen 方法始终比@Borodin 的方法慢 5 到 6 倍。最快的脚本使用了 Stream::Reader 模块。结果是一致的，每个脚本都选择了相同的前五个字符串和两个较小的文件：
100万字符核苷酸串
sequence_search_stream-reader.pl : 0.26s
sequence_search_borodin.pl : 0.39s
sequence_search_listgen.pl : 2.04s
8300万字符核苷酸串
与文件xaa中的数据：
wc xaa
0 1 83886080 xaa
% time perl sequence_search_stream-reader.pl xaa
GGCNG 31510
TAGNN 31182
AACTA 30944
GTCAN 30792
ANTAT 30756
21.33 real 20.95 user 0.35 sys
% time perl sequence_search_borodin.pl xaa
GGCNG 31510
TAGNN 31182
AACTA 30944
GTCAN 30792
ANTAT 30756
28.13 real 28.08 user 0.03 sys
% time perl sequence_search_listgen.pl xaa
GGCNG 31510
TAGNN 31182
AACTA 30944
GTCAN 30792
ANTAT 30756
157.54 real 156.93 user 0.45 sys
10亿字符核苷酸串
在更大的文件中，差异幅度相似，但是，因为它没有正确处理跨越块边界的序列，所以 List::Gen 脚本与 shell 命令行有相同的差异在此 post 的开头。较大的文件意味着许多块边界和计数差异。
sequence_search_stream-reader.pl : 252.12s
sequence_search_borodin.pl : 350.57s
sequence_search_listgen.pl : 1928.34s
块边界问题当然可以解决，但我很想知道使用 "lazy list" 方法引入的其他潜在错误或瓶颈。如果使用 slide 到 "lazily" 沿字符串移动在 CPU 用法方面有任何好处，它似乎因为需要在字符串中列出列表而变得毫无意义开始。
我对跨块边界读取数据留作实施练习并不感到惊讶（也许它无法处理 "magically"），但我想知道可能存在哪些其他 CPAN 模块或陈旧的子例程样式解决方案。
1. 在读取 TB 文件的每兆字节末尾跳过四个字符 - 因此有四个 5 字符的字符串组合意味着结果将不包括 3/10000 来自最终计数的 1%。
echo "scale=10; 100 * (1024^4/1024^2 ) * 4 / 1024^4 " | bc
.0003814697
一般来说，Perl 在逐字符处理解决方案（如上面发布的那些解决方案）方面确实很慢，但在诸如正则表达式之类的解决方案方面快得多，因为基本上你的开销主要是多少您正在执行的运算符。
因此，如果您可以将其转变为基于正则表达式的解决方案，那就更好了。
尝试这样做：
$ perl -wE 'my $str = "abcdefabcgbacbdebdbbcaebfebfebfeb"; for my $pos (0..4) { $str =~ s/^.// if $pos; say for $str =~ m/(.{5})/g }'|sort|uniq -c|sort -nr|head -n 5
3 ebfeb
2 febfe
2 bfebf
1 gbacb
1 fabcg
即我们在 $str 中有我们的字符串，然后我们传递它 5 次生成 5 个字符的序列，在第一次传递之后我们开始从字符串的前面切掉一个字符。在很多语言中，这会 真的 慢，因为你必须重新分配整个字符串，但是 perl 会针对这种特殊情况作弊，只是将字符串的索引设置为 1 + 之前是什么。
我没有对此进行基准测试，但我敢打赌，这样做比上面的算法更可行，你也可以在 perl 中进行 uniq 计数，当然是通过递增哈希（使用 / e regex 选项可能是最快的方式），但我只是将其卸载到此实现中的 |sort|uniq -c，这可能更快。
在 perl 中完成这一切的稍微改变的实现：
$ perl -wE 'my $str = "abcdefabcgbacbdebdbbcaebfebfebfeb"; my %occur; for my $pos (0..4) { substr($str, 0, 1) = "" if $pos; $occur{$_}++ for $str =~ m/(.{5})/gs }; for my $k (sort { $occur{$b} <=> $occur{$a} } keys %occur) { say "$occur{$k} $k" }'
3 ebfeb
2 bfebf
2 febfe
1 caebf
1 cgbac
1 bdbbc
1 acbde
1 efabc
1 aebfe
1 ebdbb
1 fabcg
1 bacbd
1 bcdef
1 cbdeb
1 defab
1 debdb
1 gbacb
1 bdebd
1 cdefa
1 bbcae
1 bcgba
1 bcaeb
1 abcgb
1 abcde
1 dbbca
其背后代码的漂亮格式：
my $str = "abcdefabcgbacbdebdbbcaebfebfebfeb";
my %occur;
for my $pos (0..4) {
substr($str, 0, 1) = "" if $pos;
$occur{$_}++ for $str =~ m/(.{5})/gs;
}
for my $k (sort { $occur{$b} <=> $occur{$a} } keys %occur) {
say "$occur{$k} $k";
}
即使您在处理文件之前没有将整个文件读入内存，您仍然可能 运行 内存不足。
一个 10 GiB 的文件包含将近 11E9 个序列。
如果您的序列是从一组 5 个字符中选出的 5 个字符的序列，则只有 55 = 3,125 个唯一序列，这很容易存储在内存中。
如果您的序列是从一组 5 个字符中选出的 20 个字符的序列，则有 520 = 95E12 个唯一序列，因此一个 10 GiB 文件的所有 11E9 个序列可以独一无二。那不适合记忆。
在这种情况下，我建议执行以下操作：
创建一个包含原始文件所有序列的文件。
以下是分块读取文件，而不是一次全部读取。棘手的部分是处理跨越两个块的序列。以下程序使用
sysread[1] 从文件中获取任意大小的数据块并将其附加到先前读取块的最后几个字符。最后一个细节允许对跨越块的序列进行计数。perl -e' use strict; use warnings qw( all ); use constant SEQ_LENGTH => 20; use constant CHUNK_SIZE => 1024 * 1024; my $buf = ""; while (1) { my $size = sysread(\*STDIN, $buf, CHUNK_SIZE, length($buf)); die($!) if !defined($size); last if !$size; for my $offset ( 0 .. length($buf) - SEQ_LENGTH ) { print(substr($buf, $offset, SEQ_LENGTH), "\n"); } substr($buf, 0, -(SEQ_LENGTH-1), ""); } ' <in.txt >sequences.txt对序列进行排序。
sort sequences.txt >sorted_sequences.txt计算每个序列的实例数，并将计数与序列一起存储在另一个文件中。
perl -e' use strict; use warnings qw( all ); my $last = ""; my $count; while (<>) { chomp; if ($_ eq $last) { ++$count; } else { print("$count $last\n") if $count; $last = $_; $count = 1; } } ' sorted_sequences.txt >counted_sequences.txt按计数对序列排序。
sort -rns counted_sequences.txt >sorted_counted_sequences.txt提取结果。
perl -e' use strict; use warnings qw( all ); my $last_count; while (<>) { my ($count, $seq) = split; last if $. > 5 && $count != $last_count; print("$seq $count\n"); $last_count = $count; } ' sorted_counted_sequences.txt这也打印出并列第 5 名。
这可以通过调整传递给 sort[2] 的参数来优化，但它应该提供不错的性能。
sysread比之前建议的read更快，因为后者在内部执行一系列 4 KiB 或 8 KiB 读取（取决于您的 Perl 版本）。鉴于序列的固定长度性质，您还可以将序列压缩为 ceil(log256(520)) = 6 个字节，然后将它们进行 base64 编码为 ceil(6 * 4/3) = 8 个字节。这意味着每个序列需要的字节减少 12 个，大大减少了读写量。
此答案的部分内容改编自 user:622310 根据 cc by-sa 3.0 许可的内容。