User Tools

Site Tools


n50.pl
#!/usr/bin/perl
#
# Source: 
#    http://genomics-array.blogspot.com/2011/02/calculating-n50-of-contig-assembly-file.html
#    Modified by Haeyoung Jeong
## Read Fasta File and compute length ###
#
my $length;
my $totalLength; 
my @arr;
my $num;


while(<>){
   chomp; 
   if(/^>([^\s]+)/){
       print " $length\n" unless $num == 0;
       $num++;
       print "SEQ: $1 ";
       push (@arr, $length);
       $totalLength += $length; 
       $length=0;
  } else {
  s/\s//g;
  s/\t//g;
  $length += length($_);
  }
}
push (@arr, $length); # for the last contig
print " $length\n";
$totalLength += $length;

close(FH);

my @sort = sort {$b <=> $a} @arr;
my $n50;
print "Total $num contigs ($totalLength)\n";
foreach my $val(@sort){
     $n50+=$val;
      if($n50 >= $totalLength/2){
         print "N50 contig length is $n50 and N50 value is: $val\n";
         last;
     }
}
print "Average contig length is ", $totalLength / $num, "\n";
print "Max. contig length is ", $sort[0], "\n";
n50.pl.txt · Last modified: 2021/03/17 13:09 by 127.0.0.1