view release on metacpan or  search on metacpan

lib/Bio/BPWrapper/SeqManipulations.pm  view on Meta::CPAN

sub pick_by_order {
    my ($ct, $currseq, $order_list) = @_;
    $out->write_seq($currseq) if $order_list->{$ct}
}

sub del_by_order {
    my ($ct, $currseq, $order_list) = @_; # say join "\t", @_;
    if ($order_list->{$ct}) { warn "Deleted sequence: ", $currseq->id(), "\n" }
    else { $out->write_seq($currseq) }
}

sub find_by_id {
    my ($action, $match, $currseq, $id_list) = @_;
    my $seq_id = $currseq->id();
    $filter_dispatch{$action . "_by_id"}->($match, $currseq, $id_list, $seq_id)
}

sub pick_by_id {
    my ($match, $currseq, $id_list, $seq_id) = @_;

    if ($id_list->{$seq_id}) {
        $id_list->{$seq_id}++;
        die "Multiple matches (" . $seq_id . ":" . $id_list->{$seq_id} - 1 . ") for $match found\n" if $id_list->{$seq_id} > 2;
        $out->write_seq($currseq)
    }
}

sub del_by_id {
    my ($match, $currseq, $id_list, $seq_id) = @_;

    if ($id_list->{$seq_id}) {
        $id_list->{$seq_id}++;
        warn "Deleted sequence: ", $currseq->id(), "\n"
    } else { $out->write_seq($currseq) }
}

sub find_by_re {
    my ($action, $currseq, $value) = @_;
    my $regex  = qr/$value/;
    my $seq_id = $currseq->id();
    $filter_dispatch{ $action . "_by_re" }->($currseq, $regex, $seq_id)
}

sub pick_by_re {
    my ($currseq, $regex, $seq_id) = @_;
    $out->write_seq($currseq) if $seq_id =~ /$regex/
}

sub del_by_re {
    my ($currseq, $regex, $seq_id) = @_;

    if ($seq_id =~ /$regex/) { warn "Deleted sequence: $seq_id\n" }
    else { $out->write_seq($currseq) }
}

# TODO This needs better documentation
sub find_by_ambig {
    my ($action, $currseq, $cutoff) = @_;
    my $string        = $currseq->seq();
    my $ct            = ($string =~ s/([^ATCG])/$1/gi); # won't work for AA seqs
    my $percent_ambig = $ct / $currseq->length();
    $filter_dispatch{"$action" . "_by_ambig"}->($currseq, $cutoff, $ct, $percent_ambig)
}

# TODO Probably better to change behavior when 'picking'?
sub pick_by_ambig {
    my ($currseq, $cutoff, $ct, $percent_ambig) = @_;
    $out->write_seq($currseq) if $percent_ambig > $cutoff
}

sub del_by_ambig {
    my ($currseq, $cutoff, $ct, $percent_ambig) = @_;

#    if ($percent_ambig > $cutoff) { warn "Deleted sequence: ", $currseq->id(), " number of N: ", $ct, "\n" }
    if ($ct >= $cutoff) { warn "Deleted sequence: ", $currseq->id(), " number of bad monomers: ", $ct, "\n" }
    else { $out->write_seq($currseq) }
}

sub find_by_length {
    my ($action, $currseq, $value) = @_;
    $filter_dispatch{$action . "_by_length"}->($currseq, $value)
}

sub pick_by_length {
    my ($currseq, $value) = @_;
    $out->write_seq($currseq) if $currseq->length() <= $value
}

sub del_by_length {
    my ($currseq, $value) = @_;

    if ($currseq->length() <= $value) { warn "Deleted sequence: ", $currseq->id(), " length: ", $currseq->length(), "\n" }
    else { $out->write_seq($currseq) }
}

1;
__END__

=head1 EXTENDING THIS MODULE

We encourage BioPerl developers to add command-line interface to their BioPerl methods here.

Here is how to extend.  We'll use option C<--count-codons> as an example.

=over 4

=item *
Create a new method like one of the above in the previous section.

=item *
Document your method in pod using C<=head2>. For example:

    =head2 count_codons()

    Count codons for coding sequences (e.g., a genome file consisting of
    CDS sequences). Wraps
    L<Bio::Tools::SeqStats-E<gt>count_codons()|https://metacpan.org/pod/Bio::Tools::SeqStats#count_codons>.

    =cut

See L<C<count_codons()>|/count_codons> for how this gets rendered.


=item *
Add the method to C<@EXPORT> list in C<SeqManipulations.pm>.

=item *
Add option to C<%opt_displatch> which maps the option used in C<bioaln> to the subroutine that
gets called here. For example:

   'count-codons' => \&count_codons,

=item *
Add option in to C<bioseq> script. See the code that starts: