lib/Bio/Index/Fastq.pm

   1 # BioPerl module for Bio::Index::Fastq
   2 #
   3 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   4 #
   5 # Cared for by Tony Cox <avc@sanger.ac.uk>
   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::Fastq - Interface for indexing (multiple) fastq files
  14
  15 =head1 SYNOPSIS
  16
  17     # Complete code for making an index for several
  18     # fastq files
  19     use Bio::Index::Fastq;
  20     use strict;
  21
  22     my $Index_File_Name = shift;
  23     my $inx = Bio::Index::Fastq->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 Fastq format
  30     use Bio::Index::Fastq;
  31     use strict;
  32
  33     my $Index_File_Name = shift;
  34     my $inx = Bio::Index::Fastq->new('-filename' => $Index_File_Name);
  35     my $out = Bio::SeqIO->new('-format' => 'Fastq','-fh' => \*STDOUT);
  36
  37     foreach my $id (@ARGV) {
  38         my $seq = $inx->fetch($id); # Returns Bio::Seq::Quality 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 funtionality for indexing fastq files, and
  50 retrieving the sequence from them. Note: for best results 'use strict'.
  51
  52 Bio::Index::Fastq supports the Bio::DB::BioSeqI interface, meaning
  53 it can be used as a Sequence database for other parts of bioperl
  54
  55 =head1 FEED_BACK
  56
  57 =head2 Mailing Lists
  58
  59 User feedback is an integral part of the evolution of this and other
  60 Bioperl modules. Send your comments and suggestions preferably to one
  61 of the Bioperl mailing lists.  Your participation is much appreciated.
  62
  63   bioperl-l@bioperl.org                  - General discussion
  64   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  65
  66 =head2 Support
  67
  68 Please direct usage questions or support issues to the mailing list:
  69
  70 I<bioperl-l@bioperl.org>
  71
  72 rather than to the module maintainer directly. Many experienced and
  73 reponsive experts will be able look at the problem and quickly
  74 address it. Please include a thorough description of the problem
  75 with code and data examples if at all possible.
  76
  77 =head2 Reporting Bugs
  78
  79 Report bugs to the Bioperl bug tracking system to help us keep track
  80 the bugs and their resolution.  Bug reports can be submitted via the
  81 web:
  82
  83   https://github.com/bioperl/bioperl-live/issues
  84
  85 =head1 AUTHOR - Tony Cox
  86
  87 Email - avc@sanger.ac.uk
  88
  89 =head1 APPENDIX
  90
  91 The rest of the documentation details each of the object
  92 methods. Internal methods are usually preceded with a _
  93
  94 =cut
  95
  96
  97 # Let the code begin...
  98
  99
 100 package Bio::Index::Fastq;
 101
 102 use strict;
 103
 104 use Bio::Seq;
 105
 106 use base qw(Bio::Index::AbstractSeq);
 107
 108 #
 109 # Suggested fix by Michael G Schwern <schwern@pobox.com> to
 110 # get around a clash with CPAN shell...
 111 #
 112
 113 sub _version {
 114     return 0.2;
 115 }
 116
 117 =head2 _file_format
 118
 119  Title   : _file_format
 120  Function: The file format for this package, which is needed
 121            by the SeqIO system when reading the sequence.
 122  Returns : 'Fastq'
 123
 124 =cut
 125
 126 sub _file_format {
 127     return 'Fastq';
 128 }
 129
 130
 131
 132 =head2 _index_file
 133
 134   Title   : _index_file
 135   Usage   : $index->_index_file( $file_name, $i )
 136   Function: Specialist function to index FASTQ format files.
 137             Is provided with a filename and an integer
 138             by make_index in its SUPER class.
 139   Example :
 140   Returns :
 141   Args    :
 142
 143 =cut
 144
 145 sub _index_file {
 146     my( $self,
 147         $file, # File name
 148         $i,    # Index-number of file being indexed
 149         ) = @_;
 150
 151     my( $begin,     # Offset from start of file of the start
 152                     # of the last found record.
 153         );
 154
 155     $begin = 0;
 156
 157     my $id_parser = $self->id_parser;
 158     my $c = 0;
 159     open my $FASTQ, '<', $file or $self->throw("Could not read file '$file': $!");
 160
 161     # In Windows, text files have '\r\n' as line separator, but when reading in
 162     # text mode Perl will only show the '\n'. This means that for a line "ABC\r\n",
 163     # "length $_" will report 4 although the line is 5 bytes in length.
 164     # We assume that all lines have the same line separator and only read current line.
 165     my $init_pos   = tell($FASTQ);
 166     my $curr_line  = <$FASTQ>;
 167     my $pos_diff   = tell($FASTQ) - $init_pos;
 168     my $correction = $pos_diff - length $curr_line;
 169     seek $FASTQ, $init_pos, 0; # Rewind position to proceed to read the file
 170
 171     # Main indexing loop
 172     while (<$FASTQ>) {
 173         if (/^@/) {
 174             my $begin = tell($FASTQ) - length( $_ ) - $correction;
 175             foreach my $id (&$id_parser($_)) {
 176                 $self->add_record($id, $i, $begin);
 177                 $c++;
 178             }
 179         }
 180     }
 181
 182     close $FASTQ;
 183     return ($c);
 184 }
 185
 186 =head2 id_parser
 187
 188   Title   : id_parser
 189   Usage   : $index->id_parser( CODE )
 190   Function: Stores or returns the code used by record_id to
 191             parse the ID for record from a string.  Useful
 192             for (for instance) specifying a different
 193             parser for different flavours of FASTQ file.
 194             Returns \&default_id_parser (see below) if not
 195             set. If you supply your own id_parser
 196             subroutine, then it should expect a fastq
 197             description line.  An entry will be added to
 198             the index for each string in the list returned.
 199   Example : $index->id_parser( \&my_id_parser )
 200   Returns : ref to CODE if called without arguments
 201   Args    : CODE
 202
 203 =cut
 204
 205 sub id_parser {
 206     my( $self, $code ) = @_;
 207
 208     if ($code) {
 209         $self->{'_id_parser'} = $code;
 210     }
 211     return $self->{'_id_parser'} || \&default_id_parser;
 212 }
 213
 214
 215
 216 =head2 default_id_parser
 217
 218   Title   : default_id_parser
 219   Usage   : $id = default_id_parser( $header )
 220   Function: The default Fastq ID parser for Fastq.pm
 221             Returns $1 from applying the regexp /^>\s*(\S+)/
 222             to $header.
 223   Returns : ID string
 224   Args    : a fastq header line string
 225
 226 =cut
 227
 228 sub default_id_parser {
 229     if ($_[0] =~ /^@\s*(\S+)/) {
 230         return $1;
 231     } else {
 232         return;
 233     }
 234 }
 235
 236 1;