File Coverage

blib/lib/Proch/N50.pm

Criterion	Covered	Total	%
statement	95	103	92.2
branch	19	28	67.8
condition	14	24	62.5
subroutine	11	11	100.0
pod	3	3	100.0
total	142	169	84.6

package Proch::N50;
#ABSTRACT: a small module to calculate N50 (total size, and total number of sequences) for a FASTA or FASTQ file. It's easy to install, with minimal dependencies.

use 5.012;
use warnings;
my  $opt_digits = 2;
$Proch::N50::VERSION = '1.5.7';
use File::Spec;
use JSON::PP;
use FASTX::Reader;
use File::Basename;
use Exporter qw(import);

our @EXPORT = qw(getStats getN50 jsonStats);

sub getStats {
    # Parses a FASTA/FASTQ file and returns stats
    # Parameters:
    # * filename (Str)
    # * Also return JSON string (Bool)

    my ( $file, $wantJSON, $customN ) = @_;
    if (defined $customN and ($customN > 100 or $customN < 0) ) {
      die "[Proch::N50] Custom value must be 0 < x < 100\n";
    }
    my $answer;
    $answer->{status} = 1;
    $answer->{N50}    = undef;

    # Check file existence
    # uncoverable condition right
    if ( !-e "$file" and $file ne '-' ) {
        $answer->{status}  = 0;
        $answer->{message} = "Unable to find <$file>";
    }

    # Return failed status if file not found or not readable
    if ( $answer->{status} == 0 ) {
        return $answer;
    }

    # my @aux = undef;
    my $Reader;
    if ($file ne '-') {
       $Reader = FASTX::Reader->new({ filename => "$file" });
    } else {
       $Reader = FASTX::Reader->new({ filename => '{{STDIN}}' });
    }
    my %sizes;
    my ( $n, $slen ) = ( 0, 0 );

    # Parse FASTA/FASTQ file
    while ( my $seq = $Reader->getRead() ) {
        ++$n;
        my $size = length($seq->{seq});
        $slen += $size;
        $sizes{$size}++;
    }

    my ($n50, $min, $max, $auN, $n75, $n90, $nx);
    unless ($n) {
      ($n50, $min, $max, $auN, $n75, $n90, $nx) =
      (   0,    0,    0,    0,    0,    0,   0);
       say STDERR "[n50] WARNING: Not a sequence file: $file";
    } else {
      # Invokes core _n50fromHash() routine
      ($n50, $min, $max, $auN, $n75, $n90, $nx) = _n50fromHash( \%sizes, $slen, $customN );
    }
    my $basename = basename($file);

    $answer->{N50}      = $n50 + 0;
    $answer->{N75}      = $n75 + 0;
    $answer->{N90}      = $n90 + 0;
    if (defined $customN) {
      $answer->{Ne}       = $nx  + 0;
      $answer->{"N$customN"} = $nx  + 0;
    }

    $answer->{auN}      = sprintf("%.${opt_digits}f", $auN + 0);
    $answer->{min}      = $min + 0;
    $answer->{max}      = $max + 0;
    $answer->{seqs}     = $n;
    $answer->{size}     = $slen;
    $answer->{filename} = $basename;
    $answer->{dirname}  = dirname($file);
    $answer->{path   }  = File::Spec->rel2abs(dirname($file));

    # If JSON is required return JSON
    if ( defined $wantJSON and $wantJSON ) {

        my $json = JSON::PP->new->ascii->pretty->allow_nonref;
        my $pretty_printed = $json->encode( $answer );
        $answer->{json} = $pretty_printed;

    }
    return $answer;
}

sub _n50fromHash {
    # _n50fromHash(): calculate stats from hash of lengths
    #
    # Parameters:
    # * A hash of  key={contig_length} and value={no_contigs}
    # * Sum of all contigs sizes
    my ( $hash_ref, $total_size, $custom_n ) = @_;
    my $progressive_sum = 0;
    my $auN = 0;
    my $n50 = undef;
    my $n75 = undef;
    my $n90 = undef;
    my $nx  = undef;
    my @sorted_keys = sort { $a <=> $b } keys %{$hash_ref};

    # Added in v. 0.039
    my $max =  $sorted_keys[-1];
    my $min =  $sorted_keys[0] ;
    # N50 definition: https://en.wikipedia.org/wiki/N50_statistic
    # Was '>=' in my original implementation of N50. Now complies with 'seqkit'
    # N50 Calculation

    foreach my $s ( @sorted_keys ) {
        my $ctgs_length = $s * ${$hash_ref}{$s};
        $progressive_sum +=  $ctgs_length;


        $auN += ( $ctgs_length ) * ( $ctgs_length / $total_size);

       if ( !$n50 and $progressive_sum > ( $total_size * ((100 - 50) / 100) ) ) {
         $n50 = $s;
       }

       if ( !$n75 and $progressive_sum > ( $total_size * ((100 - 75) / 100) ) ) {
         $n75 = $s;
       }
       if ( !$n90 and $progressive_sum > ( $total_size * ((100 - 90) / 100) ) ) {
         $n90 = $s;
       }
       if ( !$nx and defined $custom_n) {
         $nx = $s if ( $progressive_sum > ( $total_size * ((100 - $custom_n) / 100) ));
       }
    }
    return ($n50, $min, $max, $auN, $n75, $n90, $nx);

}

sub getN50 {

    # Invokes the full getStats returning N50 or 0 in case of error;
    my ($file) = @_;
    my $stats = getStats($file);

# Verify status and return
# uncoverable branch false
    if ( $stats->{status} ) {
        return $stats->{N50};
    } else {
        return 0;
    }
}

sub jsonStats {
  my ($file) = @_;
  my $stats = getStats($file,  'JSON');

# Return JSON object if getStats() was able to reduce one
# uncoverable branch false
  if (defined $stats->{json}) {
    return $stats->{json}
  } else {
    # Return otherwise
    return ;
  }
}


1;

__END__