File Coverage

blib/lib/BioX/Seq/Utils.pm

Criterion	Covered	Total	%
statement	43	43	100.0
branch	16	16	100.0
condition			n/a
subroutine	8	8	100.0
pod	4	4	100.0
total	71	71	100.0

package BioX::Seq::Utils;

use v5.10.1;

use strict;
use warnings;
use Exporter qw/import/;

our @EXPORT_OK = qw/
    build_ORF_regex
    all_orfs
    rev_com
    is_nucleic
/;

sub build_ORF_regex {

    my ($mode, $min_len) = @_;

    die "Missing mode"
        if (! defined $mode);
    die "Missing minimum length"
        if (! defined $min_len);

    # mode 0 : any set of codons not STOP
    # mode 1 : must end with STOP codon
    # mode 2 : must start with START codon
    # mode 3 : START -> STOP


    my $aasize = int($min_len/3);
    my $tail_size = $aasize - 1 - ($mode & 0x2 ? 1 : 0);
    my $codon = ".{3}";
    my $first_codon = $mode & 0x2 ? 'A[TU]G' : '';
    my $stop_codon = "[TU](?:AA|AG|GA|AR|RA)";
    my $last_codon  = $mode & 0x1 ? $stop_codon : '';
    my $re_orf = qr/
        \G                    # anchor at pos() - forces codon boundaries
        (?:$codon)*?          # discard codons before start codon
        (                     # begin matching ORF
        $first_codon          # match start codon
            (?! $stop_codon ) # (but only if not followed by stop codon)
        (?: $codon           # match codon(s)
            (?! $stop_codon ) # (but only if not followed by stop codon)
        ) {$tail_size,}       # any only if >= $tail_size codons long
        $codon                # add final codon
        )                     # end matching ORF
        $last_codon
    /ix;

    return $re_orf;

}

sub all_orfs {

    my ($seq, $mode, $min_len) = @_;

    my $orf_re = build_ORF_regex($mode, $min_len);

    my $len  = length $seq;
    my @orfs;
    for my $strand (0..1) {
        my $str = $strand ? rev_com($seq) : $seq;
        for my $frame (0..2) {
            pos($str) = $frame;
            while ($str =~ /$orf_re/g) {
                my ($s, $e) = map {$strand ? $len-$_+1 : $_} ($-[1]+1, $+[1]);
                push @orfs, [$1, $s, $e];
            }
        }
    }

    return @orfs;

}

sub is_nucleic {

    my ($seq) = @_;
    return $seq !~ /[^ACGTUMRWSYKVHDBN.-]/i;

}

sub rev_com {

    my ($seq) = @_;
    $seq =~ tr/Xx/Nn/;
    die "Bad input sequence" if (! is_nucleic($seq));
    $seq = reverse $seq;
    $seq =~ tr
        {ACGTMRWSYKVHDBNacgtmrwsykvhdbn}
        {TGCAKYWSRMBDHVNtgcakywsrmbdhvn};

    return $seq;

}

1;


__END__