view release on metacpan or  search on metacpan

lib/Bio/MUST/Core/SeqMask/Pmsf.pm  view on Meta::CPAN

package Bio::MUST::Core::SeqMask::Pmsf;
# ABSTRACT: Posterior mean site frequencies (PMSF) for sequence sites
$Bio::MUST::Core::SeqMask::Pmsf::VERSION = '0.252040';
use Moose;
use namespace::autoclean;

use autodie;
use feature qw(say);

use Carp;
use Const::Fast;
use List::AllUtils qw(sum each_arrayref);

extends 'Bio::MUST::Core::SeqMask';

use Bio::MUST::Core::Types;
use aliased 'Bio::MUST::Core::SeqMask::Rates';


# override superclass' Bool type
# Note: mask indices are as follow: [site][AA]
#       mask values  are freqs
has '+mask' => (
    isa => 'ArrayRef[ArrayRef[Num]]',
);

# TODO: mask non-applicable methods from superclass? (Liskov principle)

const my $PREC => 10;


sub chi_square_stats {
    my $self = shift;
    my $othr = shift;

    # check that both pmsf objects are the same length
    # potential bugs could come from constant sites etc
    my $s_width = $self->mask_len;
    my $o_width = $othr->mask_len;
    carp "[BMC] Warning: PMSF widths do not match: $s_width vs. $o_width!"
        unless $s_width == $o_width;

    my @stats;

    my $ea = each_arrayref [ $self->all_states ], [ $othr->all_states ];
    while (my ($s_freqs, $o_freqs) = $ea->() ) {
        push @stats, 0 + ( sprintf "%.${PREC}f", sum map {
            ( $o_freqs->[$_] - $s_freqs->[$_] )**2 / $s_freqs->[$_]
        } 0..$#$o_freqs );
    }   # Note: trick to get identical results across platforms
    # https://stackoverflow.com/questions/21204733/a-better-chi-square-test-for-perl

    return Rates->new( mask => \@stats );
}


# I/O methods


sub load {
    my $class  = shift;
    my $infile = shift;

    open my $in, '<', $infile;