我有一个很大的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万行还是因为我的代码而导致速度变慢,还是两者兼而有之。我正在寻找建议以加快此过程。谢谢!

最佳答案

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

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
...

关于performance - 使用嵌套的for循环的Perl脚本的性能降低,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/20850667/

10-13 01:13