view release on metacpan or  search on metacpan

bin/split-rates-ali.pl  view on Meta::CPAN

#!/usr/bin/env perl
# PODNAME: split-rates-ali.pl
# ABSTRACT: Split ALI files into subsets of sites based on site-wise statistics

use Modern::Perl '2011';
use autodie;

use Getopt::Euclid qw(:vars);
use Smart::Comments;

use Bio::MUST::Core;
use Bio::MUST::Core::Utils qw(change_suffix secure_outfile);
use aliased 'Bio::MUST::Core::Ali';
use aliased 'Bio::MUST::Core::SeqMask';


# TODO: generalize this to other executables through Utils?
my ($load, $store) = qw(load store);
my $suffix = '.ali';
my %out_args;
if ($ARGV_phylip) {
    ### Infiles (and outfiles) are in PHYLIP format
    $load  .= '_phylip';
    $store .= '_phylip';
    $suffix = '.phy';
    $out_args{chunk} = -1;
}

my $class = 'Bio::MUST::Core::SeqMask::Rates';
my $delta = 'delta_rates';
if ($ARGV_sitefreq) {
    $class = 'Bio::MUST::Core::SeqMask::Pmsf';
    $delta = 'chi_square_stats';
}

# optionally load other rates file for deltas
my $other_rates;
   $other_rates = $class->load($ARGV_other_rates) if $ARGV_other_rates;

for my $infile (@ARGV_infiles) {

    ### Processing: $infile
    my $rates = $class->load($infile);

    $infile =~ s/$_//xms for @ARGV_in_strip;
    my $alifile = change_suffix($infile, $suffix);
    my $ali = Ali->$load($alifile);
    $ali->gapify_seqs if $ARGV_from_scafos;
    ### Infile: "$alifile had " . $ali->width . ' sites'

    # optionally delete constant sites
    $ali->apply_mask( SeqMask->variable_mask($ali) ) if $ARGV_del_const;

    # optionally use delta rates instead of raw rates if available
    if ($other_rates) {
        $rates = $rates->$delta($other_rates);

        # optionally dump stats resulting from delta operation
        if ($ARGV_dump_stats) {
            my $outfile = secure_outfile(
                change_suffix($infile, '.stats'), $ARGV_out_suffix
            );
            ### Dumping site-wise stats to: $outfile
            $rates->store($outfile);
        }