perl - 计算数百 GB 数据中的子序列

Question

我正在尝试处理一个非常大的文件并计算文件中特定长度的所有序列的频率。

为了说明我在做什么，考虑一个包含序列 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用于此任务的最有效的内存工具？

Answer 1

从您自己的代码来看，您的数据文件似乎只有一行数据——没有被换行符分隔——所以我在下面的解决方案中假设了这一点。即使该行的末尾可能有一个换行符，最后五个最频繁的子序列的选择也会将其丢弃，因为它只发生一次

本程序使用 sysread 从文件中获取任意大小的数据块并将其附加到我们已经在内存中的数据中

循环体大部分与你自己的代码相似，但我使用了for的列表版本而不是 C 风格的，因为它更清晰

处理完每个chunk后，内存中的数据被截断到最后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);其中集 $chunk当我们经过while环形。这可确保正确计算跨越 2 个块的字符串。
$chunk = substr $chunk, -4语句从当前块中删除除最后四个字符之外的所有字符，以便下一次读取附加 CHUNK_SIZE从文件到剩余字符的字节数。这样搜索将继续，但从前一个块的最后 4 个字符以及下一个块开始:数据不会落入块之间的“裂缝”。