улучшить скорость парсинга fastq

@solved C# с тем же кодом в два раза быстрее

я разбираю файл phred33 fastq в perl, и это занимает значительное количество времени (порядка 15 минут). Файл fastq весит около 3 гигов. Есть ли разумные способы сделать это быстрее?

$file=shift;
open(FILE,$file);
open(FILEFA,">".$file.".fa");
open(FILEQA,">".$file.".qual");
while($line=<FILE>)
{
    chomp($line);
    if($line=~m/^@/)
    {


    $header=$line;
    $header =~ s/@/>/g;
    $seq=<FILE>;
    chomp($seq);
    $nothing=<FILE>;
    $nothing="";
    $fastq=<FILE>;

    print FILEFA $header."\n";
    print FILEFA $seq."\n";
    $seq="";
    print FILEQA $header."\n";

        @elm=split("",$fastq);
        $i=0;
        while(defined($elm[$i]))
        {
            $Q = ord($elm[$i]) - 33;
            if($Q!="-23")
            {
            print FILEQA $Q." ";
            }
            $i=$i+1;
        }
        print FILEQA "\n";
    }
}
print $file.".fa\n";
print $file.".qual\n";

person Community    schedule 24.04.2012    source источник
comment
Полагаю, если бы я был на вашем месте, я бы, наверное, запустил cachegrind...   -  person Frank Farmer    schedule 25.04.2012
comment
Небольшое гугление: biostars.org/post/show/5005/ Вы также можете попробовать прочитать с :raw : perlmonks.org/?node_id=837624 см. также stackoverflow.com/questions/1349604/ и, возможно, stackoverflow.com/questions/1052765/linux-perl-mmap-performance   -  person Frank Farmer    schedule 25.04.2012
comment
Я не очень хорошо разбираюсь в perl -- $nothing=<FILE>; читает весь файл как массив? Возможно, вы постоянно читаете весь файл?   -  person Frank Farmer    schedule 25.04.2012
comment
@FrankFarmer спасибо за ссылки, я их прочитал, но до сих пор не знаю, как исправить свой код. $nothing=‹FILE› читает одну строку :)   -  person    schedule 25.04.2012
comment
$nothing=<FILE>; $nothing = ""; содержит два ненужных присваивания. Просто используйте <FILE>.   -  person ikegami    schedule 25.04.2012


Ответы (1)


Здесь почти не используется процессор. это связано с вводом-выводом, поэтому в основном это время для чтения 3 ГБ. Есть микрооптимизация (и другие очистки), которые можно сделать.

Во-первых, всегда используйте use strict; use warnings;.

Основной код

my @elm = split(//, $fastq);
my $i=0;
while(defined($elm[$i])) {
    my $Q = ord($elm[$i]) - 33;
    if($Q!="-23") {
        print FILEQA $Q." ";
    }
    $i=$i+1;
}

Цель if($Q!="-23") состоит в том, чтобы проверить, является ли символ новой строкой, чего вам не пришлось бы делать, если бы вы сделали chomp($fastq);. (Что с кавычками вокруг -23?!)

chomp($fastq);
my @elm = split(//, $fastq);
my $i=0;
while(defined($elm[$i])) {
    my $Q = ord($elm[$i]) - 33;
    print FILEQA $Q." ";
    $i=$i+1;
}
print FILEQA "\n";

Использование цикла while только усложняет ситуацию. Используйте цикл for, если у вас есть известное количество итераций.

chomp($fastq);
for (split(//, $fastq)) {
    print FILEQA (ord($_)-33)." ";
}
print FILEQA "\n";

Это может помочь немного вывернуть это наизнанку.

chomp($fastq);
print FILEQA join(' ', map ord($_)-33, split //, $fastq), "\n";

С другой стороны, недостаточно наизнанку :)

$fastq =~ s/(.)/(ord($1)-33) . " "/eg;
print FILEQA $fastq;

Но что это мы просчитывали переводы? Тогда нам не пришлось бы повторно вызывать подпрограмму (код /e).

my %map = map { chr($_) => ($_-33)." " } 0x00..0xFF;

$fastq =~ s/(.)/$map{$1}/g;
print FILEQA $fastq;

После небольшой очистки получаем:

use strict;
use warnings;

my %map = map { chr($_) => ($_-33)." " } 0x00..0xFF;

my $file = shift;

my $fa_file   = "$file.fa";
my $qual_file = "$file.qual";

open(my $FILE,   '<', $file     ) or die $!;
open(my $FILEFA, '>', $fa_file  ) or die $!;
open(my $FILEQA, '>', $qual_file) or die $!;

while (my $header = <$FILE>) {
    next if $header !~ /^@/;

    my $seq = <$FILE>;
    <$FILE>;
    my $fastq = <$FILE>;

    $header =~ s/@/>/g;
    $fastq =~ s/(.)/$map{$1}/g;

    print $FILEFA $header;
    print $FILEFA $seq;

    print $FILEQA $header;
    print $FILEQA $fastq;
}

print "$fa_file\n";
print "$qual_file\n";
person ikegami    schedule 24.04.2012
comment
@ caseyr547, вероятно, это не будет иметь большого значения, как я объяснил. - person ikegami; 25.04.2012
comment
Возможно, вы одновременно читаете большие блоки данных (скажем, 64 КБ), используя sysread, а затем извлекаете из них строки (т.е. выполняете собственную буферизацию с большим буфером). Perl 5.14 использует считыватели 8K (и сделал их настраиваемыми во время сборки), тогда как более старый Perl использовал считыватели 4K, так что Perl 5.14 действительно поможет. - person ikegami; 25.04.2012
comment
Для любопытных: этот флаг конфигурации буфера PerlIO равен PERLIOBUF_DEFAULT_BUFSIZ. - person daxim; 25.04.2012
comment
@ikegami, я думаю, этот Perl просто медленный ... C # в два раза быстрее - person ; 25.04.2012