view release on metacpan or  search on metacpan

scripts/or.pl  view on Meta::CPAN

use Getopt::Long;
use File::Basename qw/ dirname /;
use Cwd qw/ abs_path /;
use Bio::Seq;
use Bio::SeqIO;
use Bio::ORA;

#----------------------------------------------------------
my $VERSION = '2.0';

#----------------------------------------------------------
sub hmm_disc
{
    my (
        $seq,    $name,   $id,  $translate, $evalue,
        $format, $detail, $aug, $hmm,       $filter,
        $frag,   $subset, $organism
       ) = @_;
    my $chrom;
    if    (defined $name)            { $chrom = $name . q{_} . $id; }
    elsif (defined $seq->display_id) { $chrom = $seq->display_id; }
    else                             { $chrom = 'noname_' . $id; }
    $seq->display_id($chrom);
    my $ora_obj = Bio::ORA->new($seq, $translate, $aug, $hmm);
    $ora_obj->{'_frag'} = $frag if (defined $frag);

    if (
        $ora_obj->find($evalue)
        && (
            !(defined $filter)
            || (   (defined $filter)
                && (substr($ora_obj->{'_hmmor'}[1], 2) == $filter))
           )
       )
    {

        if (defined $subset)
        {
            my $out =
              Bio::SeqIO->new('-format' => 'fasta', '-file' => ">>$subset");
            $out->write_seq($seq);
        }
        $ora_obj->show($format, $detail, $organism);
    }
    return $ora_obj->{'_verbose'};
}

sub fasta_filter
{
    my (
        $seq,    $ref,    $id,  $translate, $evalue,
        $format, $detail, $aug, $hmm,       $filter,
        $frag,   $subset, $organism
       ) = @_;
    my $mess = '* FASTA search for ' . $seq->display_id . "\n";
    my @hits = Bio::ORA->getHits($seq, 1, $ref);
    if ($#hits >= 0)
    {
        for (my $i = 0 ; $i <= $#hits ; $i++)
        {
            my ($hitstrand, $hitstart, $hitend) = split m/\|/, $hits[$i];
            my $seqstart = $hitstart - 250;
            $seqstart = 1 if ($seqstart < 1);
            my $seqend = $hitend + 250;
            $seqend = $seq->length() if ($seqend > $seq->length());
            my $seq_or = Bio::Seq->new(
                                  -seq      => $seq->subseq($seqstart, $seqend),
                                  -alphabet => 'dna',
                                  -id       => $seq->display_id . ":$seqstart-$seqend"
                                      );
            my $ora_obj = Bio::ORA->new($seq_or, $translate, $aug, $hmm);
            $ora_obj->{'_frag'} = $frag if (defined $frag);

            if (
                $ora_obj->find($evalue, $hitstrand, $seqstart, $seqend)
                && (
                    !(defined $filter)
                    || (   (defined $filter)
                        && (substr($ora_obj->{'_hmmor'}[1], 2) == $filter))
                   )
               )
            {

                if (defined $subset)
                {
                    my $out =
                      Bio::SeqIO->new('-format' => 'fasta',
                                      '-file'   => ">>$subset");
                    $out->write_seq($seq_or);
                }
                $ora_obj->show($format, $detail, $organism);
            }
            $mess .= q{ } . $ora_obj->{'_verbose'};
        }
    }
    else { $mess = '* No FASTA hit for ' . $seq->display_id . "\n"; }
    return $mess;
}

#------------------------ Main ----------------------------
my (
    $infile, $translate, $name,   $filter, $subset,
    $aug,    $frag,      $resume, $organism
   );
my ($evalue, $contigs, $verbose, $detail, $format, $out, $hmm, $ref) = (
                                   1e-10, 0, 0, 0, 'fasta', 1,
                                   abs_path(dirname(abs_path($0)) . '/or.hmm'),
                                   abs_path(dirname(abs_path($0)) . '/or.fasta')
                                                                       );
GetOptions(
           'sequence:s'  => \$infile,
           'hmmfile:s'   => \$hmm,
           'fastafile:s' => \$ref,
           'organism:s'  => \$organism,
           'c!'          => \$contigs,
           'a!'          => \$aug,
           'format:s'    => \$format,
           'expect:f'    => \$evalue,
           'name:s'      => \$name,
           'table:i'     => \$translate,
           'filter:i'    => \$filter,