使用嵌套 for 循环的 Perl 脚本性能缓慢答案

【问题标题】：Slow performance of Perl script using nested for loops使用嵌套 for 循环的 Perl 脚本性能缓慢
【发布时间】：2013-12-31 02:10:08
【问题描述】：

我有一个大的 FASTA 文件（一个基因序列，一个完整的染色体），其中每行包含 50 个字符（碱基 a、g、t 和 c）。这个文件大约有 400 万行。

我想重新组织文件，以便将每行的每个字符放在新文件的自己的行中。也就是说，将原始文件中的每 50 个字符的行转换为 50 个单字符的行。这将导致整个序列重写为单个列。最终，我希望将序列作为单列，这样我就可以放置一个相邻的列，其中包含每个碱基的基因组坐标位置。

这就是我的做法，使用 perl 并创建一组 for 循环。

unless(@ARGV) {
    # $0 name of the program being executed;
    print "\n usage: $0 filename\n\n"; 
    exit;
}

# use shift to pull off @ARGV value and return to $list;
my $fastafile = shift; 
open(FASTA, "<$fastafile");
my @count =(<FASTA>);
close FASTA;

# print scalar @count;

for ( my $i = 0; $i < scalar @count ; $i ++ ) {

#print "$count[$i]\n\n\n\n"; 
my @seq  = split( "", $count[ $i ] ); 
print " line = $i ";
for ( my $j = 0; $j < scalar @seq; $j++ ){
    #my $count =
    print "$seq[$j]  for count = $j \n"; 

    }

}

它似乎正在工作，但它很慢，非常慢。我想知道它是否因为 FASTA 文件有 400 万行而变慢，或者因为我的代码而变慢，或者两者兼而有之。我正在寻求建议以加快此过程。谢谢！

【问题讨论】：

你在用 fasta >header 做什么？
在for 循环之前，标头将被忽略。
你能解释一下你为什么要这样做吗？在我看来，这似乎是您试图以错误方式解决的问题之一，有点像您的问题是如何打开门，答案是使用钥匙，但您问的是如何从安全距离。
嗯，也许我理解错了，但输出目标实际上是“每行一个字符”还是“包含该字符的长字符串加上每行一个计数器”？
@carandraug，我最后想要的是一个双列文件，其中第一列是碱基，第二列是它的基因组坐标或碱基位置。该序列来自 UCSC 基因组浏览器。

标签： performance perl for-loop nested-loops bioinformatics

【解决方案1】：

问题在于您正在啜饮该文件。当大文件被 slurped 时，该进程将等到所有 I/O 结束后才开始处理。一个选项是逐行处理文件：

open my $fh, '<', $fastafile or die "Error opening file: $!";

while ( my $line = <$fh> ) {
    chomp $line;    # Remove the newline from the end of each line

    my @seq = split //, $line;

    # Loop from 0 to the last index of @seq
    for my $i ( 0 .. $#seq ) {
        print "$seq[$i] for count = $i\n";
    }
}

【讨论】：

你能注释( 0 .. $#seq )...这是什么意思？还是# 是错字？
关于 slurping，大约 6 秒后开始处理，因此与将输出打印到屏幕所需的时间相比，slurping 花费的时间并不短，这需要很长时间。到目前为止，我正在等待超过 30 分钟。当我的电脑风扇旋转时，终端窗口每隔一分钟左右就会暂停（死亡针轮旋转）。电脑变热了。
我明白了。我已经评论了范围部分。我无法确定在这种情况下我的代码是否会运行得更快，因为我无权访问该文件。可以试试吗？
$#var 是 Perl-ese，表示“数组 @var 中的元素数”
@BRPocock $#var 是获取@var 中最后一项的索引的语法。要获取数组中的 number 个项目，您可以使用 scalar(@var)（如果 scalar 已经在标量上下文中，则可以省略它）。

【解决方案2】：

也许以下内容会有所帮助：

use strict;
use warnings;

@ARGV or die "\n usage: $0 filename\n\n";

my $line = 0;
while (<>) {
    next if /^>/;
    chomp;

    print 'Line = ', $line++, "\n";
    my $count = 0;
    print "$_ for count = ", $count++, "\n" for split '';
    print "\n";
}

用法：perl script.pl fastaIn

上面也跳过了 fasta 标头。

样本输出：

Line = 0
T for count = 0
A for count = 1
C for count = 2
G for count = 3
A for count = 4
G for count = 5
...

【讨论】：

你能描述一下语法的含义吗？ print "$_ for count = " . $count++, "\n" for split //; 看起来很流线型。 for split // 是什么意思/做什么？
@ES55 - 不需要数组来进行你想要的处理。有问题的行首先将splits 序列转换为字符，并使用for 循环遍历每个字符——将每个字符分配给$_。 $count 在每次循环中递增。结果是一个T for count = 0，就像你在上面看到的那样。
@ES55 - 我已经将split // 更改为split ''，就像你拥有的那样。具有相同的效果。
感谢@Kenosis，我喜欢你的风格（一如既往）。尽管如此，这个过程的流程还是很慢，但我认为这是因为 FASTA 文件太大了。

【解决方案3】：

使用Bio::SeqIO 类来处理这个问题，它允许为fasta 格式设置width 和block（特定格式由Bio::SeqIO::fasta 处理）。如果我没记错的话，它有一些技巧可以处理非常大的序列，尽管我认为这些技巧仅限于写作部分（可耻的自我广告，我去年实现了其中一个）。像这样的东西应该可以正常工作：

use Bio::SeqIO;

## omit the -format option and it will try to guess the format
my $in  = Bio::SeqIO->new(-file => $fastafile, -format => 'Fasta');

while (my $seq = $in->next_seq()) {
  my $out = Bio::SeqIO->new(-file => ">outputfilename", -format => 'Fasta');
  $out->width(1); # 1 base pair per line
  $out->write_seq($seq);
}

请注意，这将如何允许在同一个文件中包含多个 fasta 序列（实验一个包含 6 个序列的 fasta 文件，用几行来感受一下）。

此外，这实际上会写入一个 real fasta 文件，因此您将无法更改代码来编写您的 2 列文件。但是你提到的问题，让第二列带有基本索引，对我来说没有多大意义。如果您知道第一个碱基的偏移量，那么第二列就是 $column_number + $offset + 1（以说明 fasta 标题）。但是 BioPerl 有办法做到这一点，请不要重新发明轮子。将序列加载为Bio::Seq 对象，并使用其方法获取子序列。

my $in  = Bio::SeqIO->new(-file => $fastafile);

while (my $seq = $in->next_seq()) {
  ## $subseq will be a string with the sequence from bp 500 to 1000
  my $subseq = $seq->subseq(500, 1000);
}

我不确定你会因此获得多少性能改进，但任何你认为可以改进的地方，请将其分享回 BioPerl 项目。

【讨论】：

【解决方案4】：

看起来你的主要限制是你打印出的数据比你读入的数据多几个数量级。

如果每行是 50 个字符 + 换行符，你“应该”写 100/51（大约两倍）的数据。

但是打印那个长字符串 "X for count = 29\n" 意味着每个输入字符要写出 15-16 个字符...

除此之外，您将消耗大量 RAM，但如今 4M 行 x 50 个字符并不是真的“太多”。不过，您不需要在这里“花费” 20M + 家务开销。

也许这是一个地方，编写自己的循环不如使用 Perl 运算符中的内置函数，如 qq aka "" ...

我还将变量构造移到循环之外，以节省构造和垃圾收集它们的时间。

 {                            # Inner scope for local $" and my vars            #"
     local $" = "\n";         # Separator character for stringifying lists      #"
     my ($line, @line);       # Avoid cons/gc during the loop
     while ($line = <$fh>)
     {
           chomp $line;       # Strip any newline
           @line = split ('', $line);
           print "@line\n";   # Stringification using $"
     }
 }

（抱歉，Stack Exchange 的语法高亮不知道 $" 是一个变量名，所以语法高亮有点奇怪。）

【讨论】：