利用perl脚本统计基因组中的GAP

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;
}
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容