上周我决定尝试Perl6并开始重新实现我的一个程序。我不得不说,Perl6对于对象编程来说非常简单,这在Perl5中对我来说非常痛苦。
我的程序必须读取和存储大文件,例如全基因组(高达3 Gb或更高,参见下面的示例1)或制表数据。
代码的第一个版本是通过逐行迭代(“genome.fa”.IO.lines)以Perl5方式制作的。对于正确的执行时间来说,它非常缓慢且无法确定。
my class fasta {
has Str $.file is required;
has %!seq;
submethod TWEAK() {
my $id;
my $s;
for $!file.IO.lines -> $line {
if $line ~~ /^\>/ {
say $id;
if $id.defined {
%!seq{$id} = sequence.new(id => $id, seq => $s);
}
my $l = $line;
$l ~~ s:g/^\>//;
$id = $l;
$s = "";
}
else {
$s ~= $line;
}
}
%!seq{$id} = sequence.new(id => $id, seq => $s);
}
}
sub MAIN()
{
my $f = fasta.new(file => "genome.fa");
}
所以在一点点RTFM之后,我改变了文件上的一个slurp,在我用for循环解析的\ n上的一个分区。这样我设法在2分钟内加载数据。好多了但还不够。作弊,我的意思是删除最多的\ n(例2),我将执行时间减少到30秒。相当不错,但并非完全满意,这种fasta格式并不是最常用的。
my class fasta {
has Str $.file is required;
has %!seq;
submethod TWEAK() {
my $id;
my $s;
say "Slurping ...";
my $f = $!file.IO.slurp;
say "Spliting file ...";
my @lines = $f.split(/\n/);
say "Parsing lines ...";
for @lines -> $line {
if $line !~~ /^\>/ {
$s ~= $line;
}
else {
say $id;
if $id.defined {
%!seq{$id} = seq.new(id => $id, seq => $s);
}
$id = $line;
$id ~~ s:g/^\>//;
$s = "";
}
}
%!seq{$id} = seq.new(id => $id, seq => $s);
}
}
sub MAIN()
{
my $f = fasta.new(file => "genome.fa");
}
所以RTFM又一次发现了语法的神奇之处。无论使用何种fasta格式,所以新版本和执行时间为45秒。不是最快的方式,而是更优雅和稳定。
my grammar fastaGrammar {
token TOP { <fasta>+ }
token fasta {<.ws><header><seq> }
token header { <sup><id>\n }
token sup { '>' }
token id { <[\d\w]>+ }
token seq { [<[ACGTNacgtn]>+\n]+ }
}
my class fastaActions {
method TOP ($/){
my @seqArray;
for $<fasta> -> $f {
@seqArray.push: seq.new(id => $f.<header><id>.made, seq => $f<seq>.made);
}
make @seqArray;
}
method fasta ($/) { make ~$/; }
method id ($/) { make ~$/; }
method seq ($/) { make $/.subst("\n", "", :g); }
}
my class fasta {
has Str $.file is required;
has %seq;
submethod TWEAK() {
say "=> Slurping ...";
my $f = $!file.IO.slurp;
say "=> Grammaring ...";
my @seqArray = fastaGrammar.parse($f, actions => fastaActions).made;
say "=> Storing data ...";
for @seqArray -> $s {
%!seq{$s.id} = $s;
}
}
}
sub MAIN()
{
my $f = fasta.new(file => "genome.fa");
}
我认为我找到了很好的解决方案来处理这些大文件,但性能仍然在Perl5之下。
作为Perl6的新手,我有兴趣知道是否有更好的方法来处理大数据,或者是否由于Perl6实现有一些限制?
作为Perl6的新手,我会问两个问题:
谢谢阅读 !
修复示例1:
>2L
CGACAATGCACGACAGAGGAAGCAGAACAGATATTTAGATTGCCTCTCATTTTCTCTCCCATATTATAGGGAGAAATATG
ATCGCGTATGCGAGAGTAGTGCCAACATATTGTGCTCTTTGATTTTTTGGCAACCCAAAATGGTGGCGGATGAACGAGAT
...
>3R
CGACAATGCACGACAGAGGAAGCAGAACAGATATTTAGATTGCCTCTCATTTTCTCTCCCATATTATAGGGAGAAATATG
ATCGCGTATGCGAGAGTAGTGCCAACATATTGTGCTCTTTGATTTTTTGGCAACCCAAAATGGTGGCGGATGAACGAGAT
...
修复示例2:
>2L
GACAATGCACGACAGAGGAAGCAGAACAGATATTTAGATTGCCTCTCAT...
>3R
TAGGGAGAAATATGATCGCGTATGCGAGAGTAGTGCCAACATATTGTGCT...
编辑我应用@Christoph和@timotimo的建议并用代码测试:
my class fasta {
has Str $.file is required;
has %!seq;
submethod TWEAK() {
say "=> Slurping / Parsing / Storing ...";
%!seq = slurp($!file, :enc<latin1>).split('>').skip(1).map: {
.head => seq.new(id => .head, seq => .skip(1).join) given .split("\n").cache;
}
}
}
sub MAIN()
{
my $f = fasta.new(file => "genome.fa");
}
该计划以2.7秒结束,非常棒!我还在小麦基因组(10 Gb)上尝试了这个代码。它完成了35.2秒。 Perl6终于不是那么慢了!
大感谢您的帮助!
一个简单的改进是使用固定宽度编码,如latin1
来加速字符解码,但我不确定这将有多大帮助。
就Rakudo的正则表达式/语法引擎而言,我发现它非常慢,所以可能确实需要采用更低级的方法。
我没有做任何基准测试,但我首先尝试的是这样的:
my %seqs = slurp('genome.fa', :enc<latin1>).split('>')[1..*].map: {
.[0] => .[1..*].join given .split("\n");
}
由于Perl6标准库是在Perl6本身中实现的,因此有时可以通过避免它来提高性能,以命令式方式编写代码,例如:
my %seqs;
my $data = slurp('genome.fa', :enc<latin1>);
my $pos = 0;
loop {
$pos = $data.index('>', $pos) // last;
my $ks = $pos + 1;
my $ke = $data.index("\n", $ks);
my $ss = $ke + 1;
my $se = $data.index('>', $ss) // $data.chars;
my @lines;
$pos = $ss;
while $pos < $se {
my $end = $data.index("\n", $pos);
@lines.push($data.substr($pos..^$end));
$pos = $end + 1
}
%seqs{$data.substr($ks..^$ke)} = @lines.join;
}
但是,如果使用的标准库的部分已经看到一些性能工作,这实际上可能会使事情变得更糟。在这种情况下,下一步将添加低级类型注释,如str
和int
,并用.index
替换对NQP builtins等例程的调用,例如nqp::index
。
如果这仍然太慢,那么你运气不好,需要切换语言,例如使用Inline::Perl5
调用Perl5或使用NativeCall
调用C.
请注意,@ timotimo已经完成了一些性能测量并写了关于它的an article。
如果我的短版本是基线,则命令式版本将性能提高2.4倍。
他实际上设法通过重写它来缩短短版本的3倍
my %seqs = slurp('genome.fa', :enc<latin-1>).split('>').skip(1).map: {
.head => .skip(1).join given .split("\n").cache;
}
最后,使用NQP内置函数重写命令式版本的速度提高了17倍,但考虑到潜在的可移植性问题,通常不鼓励编写此类代码,但如果您确实需要这种级别的性能,则可能现在需要:
use nqp;
my Mu $seqs := nqp::hash();
my str $data = slurp('genome.fa', :enc<latin1>);
my int $pos = 0;
my str @lines;
loop {
$pos = nqp::index($data, '>', $pos);
last if $pos < 0;
my int $ks = $pos + 1;
my int $ke = nqp::index($data, "\n", $ks);
my int $ss = $ke + 1;
my int $se = nqp::index($data ,'>', $ss);
if $se < 0 {
$se = nqp::chars($data);
}
$pos = $ss;
my int $end;
while $pos < $se {
$end = nqp::index($data, "\n", $pos);
nqp::push_s(@lines, nqp::substr($data, $pos, $end - $pos));
$pos = $end + 1
}
nqp::bindkey($seqs, nqp::substr($data, $ks, $ke - $ks), nqp::join("", @lines));
nqp::setelems(@lines, 0);
}