Bio/Index/Fasta.pm

   1 #
   2 #
   3 # BioPerl module for Bio::Index::Fasta
   4 #
   5 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   6 #
   7 # Cared for by James Gilbert <jgrg@sanger.ac.uk>
   8 #
   9 # You may distribute this module under the same terms as perl itself
  10
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 Bio::Index::Fasta - Interface for indexing (multiple) fasta files
  16
  17 =head1 SYNOPSIS
  18
  19     # Make an index for one or more fasta files
  20     use Bio::Index::Fasta;
  21     use strict;
  22
  23     my $Index_File_Name = shift;
  24     my $inx = Bio::Index::Fasta->new(-filename => $Index_File_Name,
  25                                      -write_flag => 1);
  26     $inx->make_index(@ARGV);
  27
  28
  29     # Once the index is made it can accessed, either in the
  30     # same script or a different one
  31     use Bio::Index::Fasta;
  32     use strict;
  33
  34     my $Index_File_Name = shift;
  35     my $inx = Bio::Index::Fasta->new(-filename => $Index_File_Name);
  36     my $out = Bio::SeqIO->new(-format => 'Fasta',
  37                               -fh => \*STDOUT);
  38
  39     foreach my $id (@ARGV) {
  40         my $seq = $inx->fetch($id); # Returns Bio::Seq object
  41          $out->write_seq($seq);
  42     }
  43
  44     # or, alternatively
  45     my $id;
  46     my $seq = $inx->get_Seq_by_id($id); # identical to fetch()
  47
  48 =head1 DESCRIPTION
  49
  50 Inherits functions for managing dbm files from Bio::Index::Abstract.pm,
  51 and provides the basic funtionallity for indexing fasta files, and
  52 retrieving the sequence from them. For best results 'use strict'.
  53
  54 Bio::Index::Fasta supports the Bio::DB::BioSeqI interface, meaning
  55 it can be used as a Sequence database for other parts of bioperl
  56
  57 Additional example code is available in scripts/index/*PLS and in
  58 the Bioperl Tutorial (L<http://www.bioperl.org/wiki/Bptutorial.pl>)
  59
  60 Note that by default the key for the sequence will be the first continuous
  61 string after the 'E<gt>' in the fasta header. If you want to use a specific
  62 substring of the fasta header you must use the id_parser() method.
  63
  64 You can also set or customize the unique key used to retrieve by
  65 writing your own function and calling the id_parser() method.
  66 For example:
  67
  68    $inx->id_parser(\&get_id);
  69    # make the index
  70    $inx->make_index($file_name);
  71
  72    # here is where the retrieval key is specified
  73    sub get_id {
  74       my $line = shift;
  75       $line =~ /^>.+gi\|(\d+)/;
  76       $1;
  77    }
  78
  79
  80 =head1 FEED_BACK
  81
  82 =head2 Mailing Lists
  83
  84 User feedback is an integral part of the evolution of this and other
  85 Bioperl modules. Send your comments and suggestions preferably to one
  86 of the Bioperl mailing lists.  Your participation is much appreciated.
  87
  88   bioperl-l@bioperl.org                  - General discussion
  89   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  90
  91 =head2 Support
  92
  93 Please direct usage questions or support issues to the mailing list:
  94
  95 I<bioperl-l@bioperl.org>
  96
  97 rather than to the module maintainer directly. Many experienced and
  98 reponsive experts will be able look at the problem and quickly
  99 address it. Please include a thorough description of the problem
 100 with code and data examples if at all possible.
 101
 102 =head2 Reporting Bugs
 103
 104 Report bugs to the Bioperl bug tracking system to help us keep track
 105 the bugs and their resolution.  Bug reports can be submitted via the
 106 web:
 107
 108   https://github.com/bioperl/bioperl-live/issues
 109
 110 =head1 AUTHOR - James Gilbert
 111
 112 Email - jgrg@sanger.ac.uk
 113
 114 =head1 APPENDIX
 115
 116 The rest of the documentation details each of the object
 117 methods. Internal methods are usually preceded with a _
 118
 119 =cut
 120
 121
 122 # Let the code begin...
 123
 124
 125 package Bio::Index::Fasta;
 126
 127 use strict;
 128 use warnings;
 129
 130 use Bio::Seq;
 131
 132 use base qw(Bio::Index::AbstractSeq);
 133
 134 #
 135 # Suggested fix by Michael G Schwern <schwern@pobox.com> to
 136 # get around a clash with CPAN shell...
 137 #
 138
 139 sub _version {
 140     return 0.2;
 141 }
 142
 143 =head2 _file_format
 144
 145  Title   : _file_format
 146  Function: The file format for this package, which is needed
 147            by the SeqIO system when reading the sequence.
 148  Returns : 'Fasta'
 149
 150 =cut
 151
 152 sub _file_format {
 153     return 'Fasta';
 154 }
 155
 156 =head2 _index_file
 157
 158   Title   : _index_file
 159   Usage   : $index->_index_file( $file_name, $i )
 160   Function: Specialist function to index FASTA format files.
 161             Is provided with a filename and an integer
 162             by make_index in its SUPER class.
 163   Example :
 164   Returns :
 165   Args    :
 166
 167 =cut
 168
 169 sub _index_file {
 170     my( $self,
 171          $file, # File name
 172          $i,    # Index-number of file being indexed
 173       ) = @_;
 174
 175     my( $begin,     # Offset from start of file of the start
 176                      # of the last found record.
 177     );
 178
 179     my $id_parser = $self->id_parser;
 180
 181     open my $FASTA, '<', $file or $self->throw("Could not read file '$file': $!");
 182
 183     # In Windows, text files have '\r\n' as line separator, but when reading in
 184     # text mode Perl will only show the '\n'. This means that for a line "ABC\r\n",
 185     # "length $_" will report 4 although the line is 5 bytes in length.
 186     # We assume that all lines have the same line separator and only read current line.
 187     my $init_pos   = tell($FASTA);
 188     my $curr_line  = <$FASTA>;
 189     my $pos_diff   = tell($FASTA) - $init_pos;
 190     my $correction = $pos_diff - length $curr_line;
 191     seek $FASTA, $init_pos, 0; # Rewind position to proceed to read the file
 192
 193     # Main indexing loop
 194     while (<$FASTA>) {
 195         if (/^>/) {
 196
 197             # the following was fixed to allow validation - cjfields
 198
 199             # $begin is the position of the first character after the '>'
 200             $begin = tell($FASTA) - length( $_ ) - $correction;
 201
 202             foreach my $id (&$id_parser($_)) {
 203                 $self->add_record($id, $i, $begin);
 204             }
 205         }
 206     }
 207     close $FASTA;
 208     return 1;
 209 }
 210
 211 =head2 id_parser
 212
 213   Title   : id_parser
 214   Usage   : $index->id_parser( CODE )
 215   Function: Stores or returns the code used by record_id to
 216             parse the ID for record from a string.  Useful
 217             for (for instance) specifying a different
 218             parser for different flavours of FASTA file.
 219             Returns \&default_id_parser (see below) if not
 220             set. If you supply your own id_parser
 221             subroutine, then it should expect a fasta
 222             description line.  An entry will be added to
 223             the index for each string in the list returned.
 224   Example : $index->id_parser( \&my_id_parser )
 225   Returns : ref to CODE if called without arguments
 226   Args    : CODE
 227
 228 =cut
 229
 230 sub id_parser {
 231     my( $self, $code ) = @_;
 232
 233     if ($code) {
 234         $self->{'_id_parser'} = $code;
 235     }
 236     return $self->{'_id_parser'} || \&default_id_parser;
 237 }
 238
 239 =head2 default_id_parser
 240
 241   Title   : default_id_parser
 242   Usage   : $id = default_id_parser( $header )
 243   Function: The default Fasta ID parser for Fasta.pm
 244             Returns $1 from applying the regexp /^>\s*(\S+)/
 245             to $header.
 246   Returns : ID string
 247   Args    : a fasta header line string
 248
 249 =cut
 250
 251 sub default_id_parser {
 252     if ($_[0] =~ /^>\s*(\S+)/) {
 253         return $1;
 254     } else {
 255         return;
 256     }
 257 }
 258
 259 1;