perl – 计算数百GB数据的子序列

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

为了说明我在做什么,请考虑一个包含序列的小输入文件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的list版本而不是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个字符开始：数据不会落入块之间的“裂缝”.

perl – 计算数百GB数据的子序列

解决方法

猜你在找的Perl相关文章