Bio/Index/Qual.pm

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