Bio/SeqIO/qual.pm

   1 # $Id$
   2 #
   3 # Copyright (c) 1997-9 bioperl, Chad Matsalla. All Rights Reserved.
   4 #           This module is free software; you can redistribute it and/or
   5 #           modify it under the same terms as Perl itself.
   6 #
   7 # Copyright Chad Matsalla
   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::SeqIO::qual - .qual file input/output stream
  16
  17 =head1 SYNOPSIS
  18
  19 Do not use this module directly.  Use it via the Bio::SeqIO class
  20 (see L<Bio::SeqIO> for details).
  21
  22   my $in_qual = Bio::SeqIO->new(-file    => $qualfile,
  23                                 -format  => 'qual',
  24                                 -width   => $width,
  25                                 -verbose => $verbose);
  26
  27 =head1 DESCRIPTION
  28
  29 This object can transform .qual (similar to fasta) objects to and from
  30 Bio::Seq::Quality objects. See L<Bio::Seq::Quality> for details.
  31
  32 Like the fasta module, it can take an argument '-width' to change the
  33 number of values per line (defaults to 50).
  34
  35 =head1 FEEDBACK
  36
  37 =head2 Mailing Lists
  38
  39 User feedback is an integral part of the evolution of this and other
  40 Bioperl modules. Send your comments and suggestions preferably to one
  41 of the Bioperl mailing lists.  Your participation is much appreciated.
  42
  43   bioperl-l@bioperl.org                  - General discussion
  44   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  45
  46 =head2 Support
  47
  48 Please direct usage questions or support issues to the mailing list:
  49
  50 L<bioperl-l@bioperl.org>
  51
  52 rather than to the module maintainer directly. Many experienced and
  53 reponsive experts will be able look at the problem and quickly
  54 address it. Please include a thorough description of the problem
  55 with code and data examples if at all possible.
  56
  57 =head2 Reporting Bugs
  58
  59 Report bugs to the Bioperl bug tracking system to help us keep track
  60 the bugs and their resolution.  Bug reports can be submitted via the web:
  61
  62   http://bugzilla.open-bio.org/
  63
  64 =head1 AUTHOR Chad Matsalla
  65
  66 Chad Matsalla
  67 bioinformatics@dieselwurks.com
  68
  69 =head1 CONTRIBUTORS
  70
  71 Jason Stajich, jason@bioperl.org
  72
  73 =head1 APPENDIX
  74
  75 The rest of the documentation details each of the object
  76 methods. Internal methods are usually preceded with a _
  77
  78 =cut
  79
  80 # Let the code begin...
  81
  82 package Bio::SeqIO::qual;
  83 use strict;
  84 use Bio::Seq::SeqFactory;
  85 use Dumpvalue;
  86
  87 my $dumper = new Dumpvalue();
  88
  89 use base qw(Bio::SeqIO);
  90
  91 our $WIDTH = 25;
  92
  93 sub _initialize {
  94     my($self,@args) = @_;
  95     $self->SUPER::_initialize(@args);
  96     my ($width) = $self->_rearrange([qw(WIDTH)], @args);
  97     $width && $self->width($width);
  98     if( ! defined $self->sequence_factory ) {
  99         $self->sequence_factory(Bio::Seq::SeqFactory->new
 100                                 (-verbose => $self->verbose(),
 101                                  -type => 'Bio::Seq::PrimaryQual'));
 102     }
 103 }
 104
 105 =head2 next_seq()
 106
 107  Title   : next_seq()
 108  Usage   : $scf = $stream->next_seq()
 109  Function: returns the next scf sequence in the stream
 110  Returns : Bio::Seq::PrimaryQual object
 111  Notes   : Get the next quality sequence from the stream.
 112
 113 =cut
 114
 115 sub next_seq {
 116     my ($self,@args) = @_;
 117     my ($qual,$seq);
 118     my $alphabet;
 119     local $/ = "\n>";
 120
 121     return unless my $entry = $self->_readline;
 122
 123     if ($entry eq '>') {        # very first one
 124         return unless $entry = $self->_readline;
 125     }
 126
 127     # original: my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s
 128     my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s
 129         or $self->throw("Can't parse entry [$entry]");
 130     my ($id,$fulldesc) = $top =~ /^\s*(\S+)\s*(.*)/
 131         or $self->throw("Can't parse fasta header");
 132     $id =~ s/^>//;
 133     # create the seq object
 134     $sequence =~ s/\n+/ /g;
 135     return $self->sequence_factory->create
 136         (-qual        => $sequence,
 137          -id         => $id,
 138          -primary_id => $id,
 139          -display_id => $id,
 140          -desc       => $fulldesc
 141         );
 142 }
 143
 144 =head2 _next_qual
 145
 146  Title   : _next_qual
 147  Usage   : $seq = $stream->_next_qual() (but do not do
 148           that. Use $stream->next_seq() instead)
 149  Function: returns the next quality in the stream
 150  Returns : Bio::Seq::PrimaryQual object
 151  Args    : NONE
 152  Notes  : An internal method. Gets the next quality in
 153                 the stream.
 154
 155 =cut
 156
 157 sub _next_qual {
 158     my $qual = next_primary_qual( $_[0], 1 );
 159     return $qual;
 160 }
 161
 162 =head2 next_primary_qual()
 163
 164  Title   : next_primary_qual()
 165  Usage   : $seq = $stream->next_primary_qual()
 166  Function: returns the next sequence in the stream
 167  Returns : Bio::PrimaryQual object
 168  Args    : NONE
 169
 170 =cut
 171
 172 sub next_primary_qual {
 173     # print("CSM next_primary_qual!\n");
 174     my( $self, $as_next_qual ) = @_;
 175     my ($qual,$seq);
 176     local $/ = "\n>";
 177
 178     return unless my $entry = $self->_readline;
 179
 180     if ($entry eq '>') {        # very first one
 181         return unless $entry = $self->_readline;
 182     }
 183
 184     my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s
 185         or $self->throw("Can't parse entry [$entry]");
 186     my ($id,$fulldesc) = $top =~ /^\s*(\S+)\s*(.*)/
 187         or $self->throw("Can't parse fasta header");
 188     $id =~ s/^>//;
 189     # create the seq object
 190     $sequence =~ s/\n+/ /g;
 191     if ($as_next_qual) {
 192         $qual = Bio::Seq::PrimaryQual->new(-qual       => $sequence,
 193                                            -id         => $id,
 194                                            -primary_id => $id,
 195                                            -display_id => $id,
 196                                            -desc       => $fulldesc
 197                                           );
 198     }
 199     return $qual;
 200 }
 201
 202
 203 =head2 width
 204
 205  Title   : width
 206  Usage   : $obj->width($newval)
 207  Function: Get/Set the number of values per line  for FASTA-like output
 208  Returns : value of width
 209  Args    : newvalue (optional)
 210
 211
 212 =cut
 213
 214 sub width{
 215    my ($self,$value) = @_;
 216    if( defined $value) {
 217       $self->{'width'} = $value;
 218     }
 219     return $self->{'width'} || $WIDTH;
 220 }
 221
 222
 223 =head2 write_seq
 224
 225  Title   : write_seq
 226  Usage   : $obj->write_seq( -source => $source,
 227                             -header  => "some information"
 228                             -oneline => 0);
 229  Function: Write out a list of quality values to a fasta-style file.
 230  Returns : Nothing.
 231  Args    : Requires a reference to a Bio::Seq::Quality object or a
 232            PrimaryQual object as the -source. Option 1: information
 233            for the header. Option 2: whether the quality score should
 234            be on a single line or not
 235  Notes   : If no -header is provided, $obj->id() will be used where
 236            $obj is a reference to either a Quality object or a
 237            PrimaryQual object. If $source->id() fails, "unknown" will
 238            be the header. If the Quality object has $source->length()
 239            of "DIFFERENT" (read the pod, luke), write_seq will use the
 240            length of the PrimaryQual object within the Quality object.
 241
 242 =cut
 243
 244 sub write_seq {
 245     my ($self,@args) = @_;
 246     my $width = $self->width;
 247     my ($source, $head, $oneline) = $self->_rearrange([qw(SOURCE HEADER ONELINE)], @args);
 248     if (!$source || ( !$source->isa('Bio::Seq::Quality') &&
 249                       !$source->isa('Bio::Seq::PrimaryQual')   )) {
 250         $self->throw("You must pass a Bio::Seq::Quality or a Bio::Seq::PrimaryQual".
 251                      " object to write_seq() as a parameter named \"source\"");
 252     }
 253     my $header = ($source->can("header") && $source->header) ?
 254         $source->header :
 255             ($source->can("id") && $source->id) ?
 256                 $source->id :
 257                     "unknown";
 258     my @quals = $source->qual();
 259     # ::dumpValue(\@quals);
 260     my $desc = $source->desc if $source->can('desc');
 261     $desc ||= '';
 262     $self->_print (">$header $desc\n");
 263     my (@slice,$max,$length);
 264     $length = $source->length();
 265
 266     if ( not(defined($oneline)) || $oneline == 0) {
 267         # $width quality values per line
 268         for (my $count = 1; $count<=$length; $count+= $width) {
 269             if ($count+$width > $length) {
 270                 $max = $length;
 271             } else {
 272                 $max = $count+$width-1;
 273             }
 274             my @slice = @{$source->subqual($count,$max)};
 275             $self->_print (join(' ',@slice), "\n");
 276         }
 277     } else {
 278         # quality values on a single line
 279         my @slice = @{$source->qual};
 280         $self->_print (join(' ',@slice), "\n");
 281     }
 282
 283     $self->flush if $self->_flush_on_write && defined $self->_fh;
 284     return 1;
 285 }
 286
 287
 288 1;
 289
 290 __END__