perl find_gap.pl genome.fa > gap.info
### find_gap.pl
local $/ = ">";
<>;
while (<>) {
s/>//;
my @a = split/\n/,$_,2;
$a[1] =~ s/\n//g;
my @seq = split//,$a[1];
foreach (0..$#seq) {
if ($seq[$_] =~ /N/) {
push @{$hash{$a[0]}},$_;
}
}
}
foreach my $chr (sort{$a cmp $b}keys %hash) {
my @gap = group(@{$hash{$chr}});
foreach (@gap) {
print "$chr\t$_->[0]\t$_->[1]\n";
}
}
sub group {
my @tmp = ();
my $tmp = $_[0];
my $cnt = 0;
my @return;
foreach (0..$#_-1) {
if ($_[$_+1] - $_[$_] != 1) {
push @{$tmp[$cnt]},$_[$_];
$cnt++;
}else{
push @{$tmp[$cnt]},$_[$_];
}
}
push @{$tmp[$cnt]},$_[-1];
foreach (@tmp) {
push @return,[(sort{$a<=>$b}@{$_})[0],scalar(@{$_})];
}
return @return;
}
