view release on metacpan or  search on metacpan

bin/prune-outliers.pl  view on Meta::CPAN

#!/usr/bin/env perl
# PODNAME: prune-outliers.pl
# ABSTRACT: Identify and discard outliers based on all-versus-all BLAST searches
# CONTRIBUTOR: Amandine BERTRAND <amandine.bertrand@doct.uliege.be>
# CONTRIBUTOR: Loic MEUNIER <loic.meunier@doct.uliege.be>

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

use File::Basename;
use File::Find::Rule;
use Getopt::Euclid qw(:vars);
use Path::Class qw(dir file);
use Smart::Comments;

use Bio::MUST::Core;
use Bio::MUST::Core::Constants qw(:dirs);
use Bio::MUST::Drivers;

use aliased 'Bio::MUST::Core::Ali';
use aliased 'Bio::MUST::Core::IdList';

my $qr_class = 'Bio::MUST::Core::Ali::Temporary';
my $db_class = 'Bio::MUST::Drivers::Blast::Database';
my $db_tmp_class = 'Bio::MUST::Drivers::Blast::Database::Temporary';

for my $indir (@ARGV_indirs) {

    ### Processing: $indir
    my @infiles = File::Find::Rule
        ->file()
        ->name( $SUFFICES_FOR{Ali} )
        ->in($indir);

    for my $infile (@infiles) {
        ### Processing: $infile

        # BLASTP/N all versus all
        my $db = $db_tmp_class->new( seqs => file($infile) );

        # Note: this is a special case where query and database are identical
        # (but autodetection won't work because query is just a filename)
        my $query = $db->filename;
        my $pgm   = $db->type eq 'nucl' ? 'blastn' : 'blastp';

        my $parser = $db->$pgm($query, {
            -evalue => $ARGV_evalue,
            -outfmt => 6,
            $db->type eq 'nucl' ? ( -task => 'blastn' ) : ()
        } );

        # parsing
        my %count_for;
        while ( my $hit = $parser->next_hit ) {
            my $query_id = $db->long_id_for( $hit->query_id );
            $count_for{$query_id}++;
        }

        my $outdir = dir($indir)->basename . '-pruned';

        for (my $t = $ARGV_min_ident; $t <= $ARGV_max_ident; $t += 0.1) {

            my @ids
                = grep { ($count_for{$_} / keys %count_for) >= $t }
                  grep {  $count_for{$_} >= $ARGV_min_hits } keys %count_for
            ;