File Coverage

blib/lib/Bio/Perl.pm

Criterion	Covered	Total	%
statement	101	179	56.4
branch	26	82	31.7
condition	2	14	14.2
subroutine	16	21	76.1
pod	13	13	100.0
total	158	309	51.1

package Bio::Perl;
$Bio::Perl::VERSION = '1.7.4';
use vars qw(@EXPORT @EXPORT_OK $DBOKAY);
use strict;
use Carp;

use Bio::SeqIO;
use Bio::Seq;
use Bio::Root::Version '$VERSION';
BEGIN {
    eval {
        require Bio::DB::EMBL;
        require Bio::DB::GenBank;
        require Bio::DB::SwissProt;
        require Bio::DB::RefSeq;
        require Bio::DB::GenPept;
    };
    if( $@ ) {
        $DBOKAY = 0;
    } else {
        $DBOKAY = 1;
    }
}

use base qw(Exporter);

@EXPORT = qw(read_sequence read_all_sequences write_sequence
             new_sequence get_sequence translate translate_as_string
             reverse_complement revcom revcom_as_string
             reverse_complement_as_string blast_sequence write_blast);

@EXPORT_OK = @EXPORT;



sub read_sequence{
    my ($filename,$format) = @_;

    if( !defined $filename ) {
        confess "read_sequence($filename) - usage incorrect";
    }

    my $seqio;

    if( defined $format ) {
        $seqio = Bio::SeqIO->new( '-file' => $filename, '-format' => $format);
    } else {
        $seqio = Bio::SeqIO->new( '-file' => $filename);
    }

    my $seq = $seqio->next_seq();

    return $seq;
}



sub read_all_sequences{
    my ($filename,$format) = @_;

    if( !defined $filename ) {
        confess "read_all_sequences($filename) - usage incorrect";
    }

    my $seqio;

    if( defined $format ) {
        $seqio = Bio::SeqIO->new( '-file' => $filename, '-format' => $format);
    } else {
        $seqio = Bio::SeqIO->new( '-file' => $filename);
    }

    my @seq_array;

    while( my $seq = $seqio->next_seq() ) {
        push(@seq_array,$seq);
    }

    return @seq_array;
}



sub write_sequence{
    my ($filename,$format,@sequence_objects) = @_;

    if( scalar(@sequence_objects) == 0 ) {
        confess("write_sequence(filename,format,sequence_object)");
    }

    my $error = 0;
    my $seqname = "sequence1";

    # catch users who haven't passed us a filename we can open
    if( $filename !~ /^\>/ && $filename !~ /^|/ ) {
        $filename = ">".$filename;
    }

    my $seqio = Bio::SeqIO->new('-file' => $filename, '-format' => $format);

    foreach my $seq ( @sequence_objects ) {
        my $seq_obj;

        if( !ref $seq ) {
            if( length $seq > 50 ) {
                # odds are this is a sequence as a string, and someone has not figured out
                # how to make objects. Warn him/her and then make a sequence object
                # from this
                if( $error == 0 ) {
                    carp("WARNING: You have put in a long string into write_sequence.\n".
                         "I suspect this means that this is the actual sequence\n".
                         "In the future try the\n".
                         "  new_sequence method of this module to make a new sequence object.\n".
                         "Doing this for you here\n");
                    $error = 1;
                }

                $seq_obj = new_sequence($seq,$seqname);
                $seqname++;
            } else {
                confess("You have a non object [$seq] passed to write_sequence. It maybe that you".
                        "want to use new_sequence to make this string into a sequence object?");
            }
        } else {
            if( !$seq->isa("Bio::SeqI") ) {
                confess("object [$seq] is not a Bio::Seq object; can't write it out");
            }
            $seq_obj = $seq;
        }

        # finally... we get to write out the sequence!
        $seqio->write_seq($seq_obj);
    }
    1;
}


sub new_sequence{
    my ($seq,$name,$accession) = @_;

    if( !defined $seq ) {
        confess("new_sequence(sequence_as_string) usage");
    }

    $name ||= "no-name-for-sequence";

    my $seq_object = Bio::Seq->new( -seq => $seq, -id => $name);

    $accession && $seq_object->accession_number($accession);

    return $seq_object;
}


sub blast_sequence {
    my ($seq,$verbose) = @_;

    if( !defined $verbose ) {
        $verbose = 1;
    }

    if( !ref $seq ) {
        $seq = Bio::Seq->new( -seq => $seq, -id => 'blast-sequence-temp-id');
    } elsif ( !$seq->isa('Bio::PrimarySeqI') ) {
        croak("[$seq] is an object, but not a Bio::Seq object, cannot be blasted");
    }

    require Bio::Tools::Run::RemoteBlast;

    my $prog = ( $seq->alphabet eq 'protein' ) ? 'blastp' : 'blastn';
    my $e_val= '1e-10';

    my @params = ( '-prog' => $prog,
                   '-expect' => $e_val,
                   '-readmethod' => 'SearchIO' );

    my $factory = Bio::Tools::Run::RemoteBlast->new(@params);

    my $r = $factory->submit_blast($seq);
    if( $verbose ) {
        print STDERR "Submitted Blast for [".$seq->id."] ";
    }
    sleep 5;

    my $result;

    LOOP :
    while( my @rids = $factory->each_rid) {
        foreach my $rid ( @rids ) {
            my $rc = $factory->retrieve_blast($rid);
            if( !ref($rc) ) {
                if( $rc < 0 ) {
                    $factory->remove_rid($rid);
                }
                if( $verbose ) {
                    print STDERR ".";
                }
                sleep 10;
            } else {
                $result = $rc->next_result();
                $factory->remove_rid($rid);
                last LOOP;
            }
        }
    }

    if( $verbose ) {
        print STDERR "\n";
    }
    return $result;
}


sub write_blast {
    my ($filename,$blast) = @_;

    if( $filename !~ /^\>/ && $filename !~ /^|/ ) {
        $filename = ">".$filename;
    }

    my $output = Bio::SearchIO->new( -output_format => 'blast', -file => $filename);

    $output->write_result($blast);

}


my $genbank_db = undef;
my $genpept_db = undef;
my $embl_db = undef;
my $swiss_db = undef;
my $refseq_db = undef;

sub get_sequence{
    my ($db_type,$identifier) = @_;
    if( ! $DBOKAY ) {
        confess ("Your system does not have one of LWP, HTTP::Request::Common, IO::String\n".
                 "installed so the DB retrieval method is not available.\n".
                 "Full error message is:\n $!\n");
        return;
    }
    $db_type = lc($db_type);

    my $db;

    if( $db_type =~ /genbank/ ) {
        if( !defined $genbank_db ) {
            $genbank_db = Bio::DB::GenBank->new();
        }
        $db = $genbank_db;
    }
    if( $db_type =~ /genpept/ ) {
        if( !defined $genpept_db ) {
            $genpept_db = Bio::DB::GenPept->new();
        }
        $db = $genpept_db;
    }

    if( $db_type =~ /swiss/ ) {
        if( !defined $swiss_db ) {
            $swiss_db = Bio::DB::SwissProt->new();
        }
        $db = $swiss_db;
    }

    if( $db_type =~ /embl/ ) {
        if( !defined $embl_db ) {
            $embl_db = Bio::DB::EMBL->new();
        }
        $db = $embl_db;
    }

    if( $db_type =~ /refseq/ or ($db_type !~ /swiss/ and
                                 $identifier =~ /^\s*N\S+_/)) {
        if( !defined $refseq_db ) {
            $refseq_db = Bio::DB::RefSeq->new();
        }
        $db = $refseq_db;
    }

    my $seq;

    if( $identifier =~ /^\w+\d+$/ ) {
        $seq = $db->get_Seq_by_acc($identifier);
    } else {
        $seq = $db->get_Seq_by_id($identifier);
    }

    return $seq;
}



sub translate {
    my ($scalar) = shift;

    my $obj;
    if( ref $scalar ) {
        if( !$scalar->isa("Bio::PrimarySeqI") ) {
            confess("Expecting a sequence object not a $scalar");
        } else {
            $obj= $scalar;
        }
    } else {
        # check this looks vaguely like DNA
        my $n = ( $scalar =~ tr/ATGCNatgcn/ATGCNatgcn/ );
        if( $n < length($scalar) * 0.85 ) {
            confess("Sequence [$scalar] is less than 85% ATGCN, which doesn't look very DNA to me");
        }
        $obj = Bio::PrimarySeq->new(-id => 'internalbioperlseq',-seq => $scalar);
    }
    return $obj->translate();
}



sub translate_as_string {
    my ($scalar) = shift;
    my $obj = Bio::Perl::translate($scalar);
    return $obj->seq;
}



sub reverse_complement {
    my ($scalar) = shift;

    my $obj;

    if( ref $scalar ) {
        if( !$scalar->isa("Bio::PrimarySeqI") ) {
            confess("Expecting a sequence object not a $scalar");
        } else {
            $obj= $scalar;
        }

    } else {

        # check this looks vaguely like DNA
        my $n = ( $scalar =~ tr/ATGCNatgcn/ATGCNatgcn/ );

        if( $n < length($scalar) * 0.85 ) {
            confess("Sequence [$scalar] is less than 85% ATGCN, which doesn't look very DNA to me");
        }

        $obj = Bio::PrimarySeq->new(-id => 'internalbioperlseq',-seq => $scalar);
    }

    return $obj->revcom();
}


sub revcom {
    return &Bio::Perl::reverse_complement(@_);
}



sub reverse_complement_as_string {
    my ($scalar) = shift;
    my $obj = &Bio::Perl::reverse_complement($scalar);
    return $obj->seq;
}



sub revcom_as_string {
    my ($scalar) = shift;
    my $obj = &Bio::Perl::reverse_complement($scalar);
    return $obj->seq;
}


1;


# ABSTRACT: Functional access to BioPerl for people who don't know objects
# AUTHOR: Ewan Birney
# OWNER: many people (see the individual modules for their copyright holders)
# LICENSE: Perl_5

__END__