我有一个很大的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/