############################################################################# # $Author: markus $ # $Date: 2009-11-12 12:42:02 +0100 (Thu, 12 Nov 2009) $ # $Revision: 1842 $ ############################################################################# package Bio::Grep; use strict; use warnings; require UNIVERSAL::require; use base 'Bio::Root::Root'; use version; our $VERSION = qv('0.10.6'); sub new { my ( $class, $backendname ) = @_; my $self = $class->SUPER::new(); if ( !defined $backendname ) { $backendname = 'Vmatch'; } my %known_backends = ( Agrep => 1, Vmatch => 1, GUUGle => 1, RE => 1, ); if ( !defined $known_backends{$backendname} ) { $self->throw( -class => 'Bio::Root::BadParameter', -text => 'Unknown back-end.', -value => $backendname . ' not supported.' ); } my $module = "Bio::Grep::Backend::$backendname"; $module->require; my $backend = $module->new(); $backend->settings->tmppath( File::Spec->tmpdir() ); return $backend; } 1; # Magic true value required at end of module __END__ =head1 NAME Bio::Grep - Perl extension for searching in DNA and Protein sequences =head1 VERSION This document describes Bio::Grep version 0.10.6 =head1 SYNOPSIS use Bio::Grep; my $sbe = Bio::Grep->new('Vmatch'); # define the location of the suffix arrays $sbe->settings->datapath('data'); mkdir($sbe->settings->datapath); # now generate a suffix array. you have to do this only once. $sbe->generate_database({ file => 'ATH1.cdna', description => 'AGI Transcripts', }); # search in this suffix array $sbe->settings->database('ATH1.cdna'); # search for the reverse complement and allow 2 mismatches $sbe->settings->query('UGAACAGAAAG'); $sbe->settings->reverse_complement(1); $sbe->settings->mismatches(2); # or you can use Fasta file with queries # $sbe->settings->query_file('Oligos.fasta'); # $sbe->search(); # Alternatively, you can specify the settings in the search call. # This also resets everything except the paths and the database # (because it is likely that they don't change when search is called # multiple times) $sbe->search( { query => 'AGAGCCCT', reverse_complement => 1, mismatches => 1, }); my @ids; while ( my $res = $sbe->next_res ) { print $res->sequence->id . "\n"; print $res->alignment_string() . "\n\n"; push @ids, $res->sequence_id; } # get the gene sequences of all matches as Bio::SeqIO object. # (to generate a Fasta file for example) my $seqio = $sbe->get_sequences(\@ids); =head1 DESCRIPTION B is a collection of Perl modules for searching in DNA and Protein sequences. It supports different back-ends, most importantly some (enhanced) suffix array implementations. Currently, there is no suffix array tool that works in all scenarios (for example whole genome, protein and RNA data). C provides a common API to the most popular tools. This way, you can easily switch or combine tools. =head1 METHODS =head2 CONSTRUCTOR =over =item C This method constructs a C back-end object. Available external back-ends are C, C, and C. Perl regular expressions are available in the C back-end. C is default. Sets temporary path to Ctmpdir();> my $sbe = Bio::Grep->new('Agrep'); Returns an object that uses L as base class. See L, L, L, L and L. =back =head1 FEATURES =over =item C supports most of the features of the back-ends. If you need a particular feature that is not supported, write a feature request. In general it should be easy to integrate. For a complete list of supported features, see L, for an overview see L<"FEATURE COMPARISON">. =item This module should be suitable for large data sets. The back-end output is piped to a temporary file and the parser only stores the current hit in memory. =item C includes an interface for search result filters. See L<"FILTERS">. This module also allows you to retrieve up- and downstream regions. Together with filters, this makes C an ideal framework for seed and extend algorithms. =item C was in particular designed for web services and therefore checks the settings carefully before calling back-ends. See L<"SECURITY">. =back =head1 QUICK START This is only a short overview of the functionality of this module. You should also read L and the documentation of the back-end you want to use (e.g. L). L is a (not yet comprehensive) collection of recipes for common problems. =head2 GENERATE DATABASES As a first step, you have to generate a C database out of your Fasta file in which you want to search. A C database consists of a couple of files and allows you to retrieve information about the database as well as to perform queries as fast and memory efficient as possible. You have to do this only once for every file. For example: my $sbe = Bio::Grep->new('Vmatch'); $sbe->generate_database({ file => 'ATH1.cdna', datapath => 'data', description => 'AGI Transcripts', }); Now, in a second script: my $sbe = Bio::Grep->new('Vmatch'); $sbe->settings->datapath('data'); my %local_dbs_description = $sbe->get_databases(); my @local_dbs = sort keys %local_dbs_description; Alternatively, you can use L which is part of this distribution: bgrep --backend Vmatch --database TAIR6_cdna_20060907 --datapath data --createdb =head2 SEARCH SETTINGS All search settings are stored in the L object of the back-end: $sbe->settings To set an option, call $sbe->settings->optionname(value) For example $sbe->settings->datapath('data'); # take first available database $sbe->settings->database($local_dbs[0]); $sbe->settings->query('AGAGCCCT'); See the documentation of your back-end for available options. =head2 SEARCH To start the back-end with the specified settings, simply call $sbe->search(); This method also accepts an hash reference with settings. In this case, all previous defined options except all paths and the database are set to their default values. $sbe->search({ mismatches => 2, reverse_complement => 0, query => 'AGAGCCCT' }); =head2 ANALYZE SEARCH RESULTS Use such a L like while loop to analyze the search results. while ( my $res = $sbe->next_res ) { print $res->sequence->id . "\n"; print $res->alignment_string() . "\n\n"; } See L for all available information. =head1 BGREP This distribution comes with a sample script called L. =head1 WHICH BACK-END? We support these external back-ends: =over =item C L =item C L (original Wu-Manber 1992 implementation for UNIX), L (DOS, Windows, OS/2), L (Agrep binary of C) and L (TRE implementation). =item C L =back =head2 FEATURE COMPARISON =begin html

Feature	Agrep	GUUGle	RE	Vmatch
Suffix Arrays/Trees	no	yes	no	yes
Sliding Window	yes	no	yes	no
Persistent Index¹	no	no	no	yes
Mismatches	yes	no	no	yes
Edit Distance	yes	no	no	yes
Insertions	no	no	no	no
Deletions	no	no	no	no
Multiple Queries²	no	yes	no	yes
GU	no	yes	no	no
DNA/RNA	yes	yes	yes	yes
Protein	yes	no	yes	yes
Direct and Revcom	no	yes	yes	yes
Reverse Complement	yes	yes	yes	yes
Upstream/Downstream	no	yes	yes	yes
Filters	no	yes	yes	yes
Query Length³	no	yes	no	yes
Regular Expressions⁴	no	no	yes	no

¹Needs pre-calculation and (much) more memory but queries are in general faster
²With query_file
³Matches if a substring of the query of size n or larger matches
⁴Agrep soon

=end html =begin man Features || Agrep | GUUGle | RE | Vmatch Suffix Arrays/Trees || no | yes | no | yes Sliding Window || yes | no | yes | no Persistent Index 1 || no | no | no | yes Mismatches || yes | no | no | yes Edit Distance || yes | no | no | yes Insertions || no | no | no | no Deletions || no | no | no | no Multiple Queries 2 || no | yes | no | yes GU || no | yes | no | no DNA/RNA || yes | yes | yes | yes Protein || yes | no | yes | yes Direct and Revcom || no | yes | yes | yes Reverse Complement || yes | yes | yes | yes Upstream/Downstream || no | yes | yes | yes Filters || no | yes | yes | yes Query Length 3 || no | yes | no | yes Regular Expressions 4 || no | no | yes | no -- 1 Needs pre-calculation and (much) more memory but queries are in general faster 2 With query_file 3 Matches if a substring of the query of size n or larger matches 4 Agrep soon =end man C is fast but needs a lot of memory. C is the best choice if you allow many mismatches in short sequences, if you want to search in Fasta files with relatively short sequences (e.g CDNA or Protein databases) and if you are only interested in which sequences the approximate match was found. Its performance is in this case amazing. If you want the exact positions of a match in the sequence, choose C. If you want nice alignments, choose C too (C can automatically align the sequence and the query in the C back-end, but then C is faster). Filters require exact positions, so you can't use them with C. This may change in future version or not. The C implementation of the C library (L) is also supported. This implementation has less limitations and more features (e.g. you get the exact hit positions) but is much slower. See L. C may be the best choice if you have RNA queries (counts GU as no mismatch) and if you are interested in only exact matches. Another solution here would be to use C and write a filter (see next section) that only allows GU mismatches. Of course, this is only an alternative if you can limit (C<$sbe-Esettings-Emismatches()>) the maximal number of GU mismatches. C with its pre-calculated suffix arrays is really fast, so you should consider this option. Perl regular expressions are available in the C back-end. It is a very simple back-end which is written in pure Perl and which does not require any additional software. =head1 FILTERS Filtering search results is a common task. For that, C provides an filter interface, L. Writing filters is straightforward: package MyFilter; use strict; use warnings; use Bio::Grep::Filter::FilterI; use base 'Bio::Grep::Filter::FilterI'; use Class::MethodMaker [ new => [ qw / new2 / ], ... # here some local variables, see perldoc Class::MethodMaker ]; sub new { my $self = shift->new2; $self->delete(1); # a filter that actually filters, not only adds # remarks to $self->search_result->remark $self->supports_alphabet( dna => 1, protein => 1); $self; } sub filter { my $self = shift; # code that examines $self->search_result # and returns 0 (not passed) or 1 (passed) ... $self->message('passed'); return 1; } sub reset_filter { my $self = shift; # if you need local variables, you can clean up here } 1;# Magic true value required at end of module To apply your filter: ... my $filter = MyFilter->new(); $sbe->settings->filters( ( $filter ) ); $sbe->search(); See L. =head1 ERROR HANDLING C throws L exceptions when errors occur. You can use the module L to catch these exceptions: use Error qw(:try); ... try { $sbe->search(); } catch Bio::Root::SystemException with { my $E = shift; print STDERR 'Back-end call failed: ' . $E->{'-text'} . ' (' . $E->{'-line'} . ")\n"; exit(1); } catch Bio::Root::BadParameter with { my $E = shift; print STDERR 'Wrong Settings: ' . $E->{'-text'} . ' (' . $E->{'-line'} . ")\n"; exit(1); } otherwise { my $E = shift; print STDERR "An unexpected exception occurred: \n$E"; exit(1); }; C throws a C when a system() call failed, C whenever C recognizes some problems in the settings. Be aware that C does not find all of these problems. In such a case, the back-end call will fail and this module will throw a C. Whenever it is not possible to open, copy, close, delete or write a file, croak() (L) is called. See L, L. =head1 SECURITY The use of this module (in Web Services for example) should be quite secure. All test run in taint mode. C checks the settings before it generates the string for the system() call and uses L for all temporary files. However, keep in mind that it is quite B without any further settings check, especially with the C back-end. So you should limit C, C and all other settings that have an significant impact on the calculation time. You should also set C. =head1 INCOMPATIBILITIES None reported. =head1 BUGS AND LIMITATIONS No bugs have been reported. There is not yet a nice interface for searching for multiple queries. However, C and C support this feature. So you can generate a Fasta query file with L and then set C<$sbe-Esettings-Equery_file()>. To find out, to which query a match belongs, you have to check C<$res-Equery>. It is likely that C<$sbe-Esettings-Equery> is renamed to C. Please report any bugs or feature requests to C, or through the web interface at L. =head1 SEE ALSO L L L L L L L =head2 PUBLICATIONS C: L =head1 AUTHOR Markus Riester, Emriester@gmx.deE =head1 LICENSE AND COPYRIGHT Copyright (C) 2007-2009 by M. Riester. Based on Weigel::Search v0.13, Copyright (C) 2005-2006 by Max Planck Institute for Developmental Biology, Tuebingen. This module is free software; you can redistribute it and/or modify it under the same terms as Perl itself. =head1 DISCLAIMER OF WARRANTY BECAUSE THIS SOFTWARE IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY FOR THE SOFTWARE, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE SOFTWARE "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE SOFTWARE IS WITH YOU. SHOULD THE SOFTWARE PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR, OR CORRECTION. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE THE SOFTWARE AS PERMITTED BY THE ABOVE LICENSE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE SOFTWARE (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE SOFTWARE TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. =cut # vim: ft=perl sw=4 ts=4 expandtab