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