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

为了说明我在做什么,请考虑一个包含序列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);在我们通过$chunk循环时设置while。这样可以确保跨2个块的字符串正确计数。

$chunk = substr $chunk, -4语句从当前块中删除除最后四个字符以外的所有字符,以便下一次读取将文件中的CHUNK_SIZE字节追加到其余字符。这样,搜索将继续,但从下一个块之外的前一个块的字符的后4个字符开始:数据不会落入块之间的“裂缝”中。

07-24 19:59
查看更多