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 Reporting Bugs
  47
  48 Report bugs to the Bioperl bug tracking system to help us keep track
  49 the bugs and their resolution.  Bug reports can be submitted via the web:
  50
  51   http://bugzilla.open-bio.org/
  52
  53 =head1 AUTHOR Chad Matsalla
  54
  55 Chad Matsalla
  56 bioinformatics@dieselwurks.com
  57
  58 =head1 CONTRIBUTORS
  59
  60 Jason Stajich, jason@bioperl.org
  61
  62 =head1 APPENDIX
  63
  64 The rest of the documentation details each of the object
  65 methods. Internal methods are usually preceded with a _
  66
  67 =cut
  68
  69 # Let the code begin...
  70
  71 package Bio::SeqIO::qual;
  72 use strict;
  73 use Bio::Seq::SeqFactory;
  74 use Dumpvalue;
  75
  76 my $dumper = new Dumpvalue();
  77
  78 use base qw(Bio::SeqIO);
  79
  80 our $WIDTH = 25;
  81
  82 sub _initialize {
  83     my($self,@args) = @_;
  84     $self->SUPER::_initialize(@args);
  85     my ($width) = $self->_rearrange([qw(WIDTH)], @args);
  86     $width && $self->width($width);
  87     if( ! defined $self->sequence_factory ) {
  88         $self->sequence_factory(Bio::Seq::SeqFactory->new
  89                                 (-verbose => $self->verbose(),
  90                                  -type => 'Bio::Seq::PrimaryQual'));
  91     }
  92 }
  93
  94 =head2 next_seq()
  95
  96  Title   : next_seq()
  97  Usage   : $scf = $stream->next_seq()
  98  Function: returns the next scf sequence in the stream
  99  Returns : Bio::Seq::PrimaryQual object
 100  Notes   : Get the next quality sequence from the stream.
 101
 102 =cut
 103
 104 sub next_seq {
 105     my ($self,@args) = @_;
 106     my ($qual,$seq);
 107     my $alphabet;
 108     local $/ = "\n>";
 109
 110     return unless my $entry = $self->_readline;
 111
 112     if ($entry eq '>') {        # very first one
 113         return unless $entry = $self->_readline;
 114     }
 115
 116     # original: my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s
 117     my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s
 118         or $self->throw("Can't parse entry [$entry]");
 119     my ($id,$fulldesc) = $top =~ /^\s*(\S+)\s*(.*)/
 120         or $self->throw("Can't parse fasta header");
 121     $id =~ s/^>//;
 122     # create the seq object
 123     $sequence =~ s/\n+/ /g;
 124     return $self->sequence_factory->create
 125         (-qual        => $sequence,
 126          -id         => $id,
 127          -primary_id => $id,
 128          -display_id => $id,
 129          -desc       => $fulldesc
 130         );
 131 }
 132
 133 =head2 _next_qual
 134
 135  Title   : _next_qual
 136  Usage   : $seq = $stream->_next_qual() (but do not do
 137           that. Use $stream->next_seq() instead)
 138  Function: returns the next quality in the stream
 139  Returns : Bio::Seq::PrimaryQual object
 140  Args    : NONE
 141  Notes  : An internal method. Gets the next quality in
 142                 the stream.
 143
 144 =cut
 145
 146 sub _next_qual {
 147     my $qual = next_primary_qual( $_[0], 1 );
 148     return $qual;
 149 }
 150
 151 =head2 next_primary_qual()
 152
 153  Title   : next_primary_qual()
 154  Usage   : $seq = $stream->next_primary_qual()
 155  Function: returns the next sequence in the stream
 156  Returns : Bio::PrimaryQual object
 157  Args    : NONE
 158
 159 =cut
 160
 161 sub next_primary_qual {
 162     # print("CSM next_primary_qual!\n");
 163     my( $self, $as_next_qual ) = @_;
 164     my ($qual,$seq);
 165     local $/ = "\n>";
 166
 167     return unless my $entry = $self->_readline;
 168
 169     if ($entry eq '>') {        # very first one
 170         return unless $entry = $self->_readline;
 171     }
 172
 173     my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s
 174         or $self->throw("Can't parse entry [$entry]");
 175     my ($id,$fulldesc) = $top =~ /^\s*(\S+)\s*(.*)/
 176         or $self->throw("Can't parse fasta header");
 177     $id =~ s/^>//;
 178     # create the seq object
 179     $sequence =~ s/\n+/ /g;
 180     if ($as_next_qual) {
 181         $qual = Bio::Seq::PrimaryQual->new(-qual       => $sequence,
 182                                            -id         => $id,
 183                                            -primary_id => $id,
 184                                            -display_id => $id,
 185                                            -desc       => $fulldesc
 186                                           );
 187     }
 188     return $qual;
 189 }
 190
 191
 192 =head2 width
 193
 194  Title   : width
 195  Usage   : $obj->width($newval)
 196  Function: Get/Set the number of values per line  for FASTA-like output
 197  Returns : value of width
 198  Args    : newvalue (optional)
 199
 200
 201 =cut
 202
 203 sub width{
 204    my ($self,$value) = @_;
 205    if( defined $value) {
 206       $self->{'width'} = $value;
 207     }
 208     return $self->{'width'} || $WIDTH;
 209 }
 210
 211
 212 =head2 write_seq
 213
 214  Title   : write_seq
 215  Usage   : $obj->write_seq( -source => $source,
 216                             -header  => "some information"
 217                             -oneline => 0);
 218  Function: Write out a list of quality values to a fasta-style file.
 219  Returns : Nothing.
 220  Args    : Requires a reference to a Bio::Seq::Quality object or a
 221            PrimaryQual object as the -source. Option 1: information
 222            for the header. Option 2: whether the quality score should
 223            be on a single line or not
 224  Notes   : If no -header is provided, $obj->id() will be used where
 225            $obj is a reference to either a Quality object or a
 226            PrimaryQual object. If $source->id() fails, "unknown" will
 227            be the header. If the Quality object has $source->length()
 228            of "DIFFERENT" (read the pod, luke), write_seq will use the
 229            length of the PrimaryQual object within the Quality object.
 230
 231 =cut
 232
 233 sub write_seq {
 234     my ($self,@args) = @_;
 235     my $width = $self->width;
 236     my ($source, $head, $oneline) = $self->_rearrange([qw(SOURCE HEADER ONELINE)], @args);
 237     if (!$source || ( !$source->isa('Bio::Seq::Quality') &&
 238                       !$source->isa('Bio::Seq::PrimaryQual')   )) {
 239         $self->throw("You must pass a Bio::Seq::Quality or a Bio::Seq::PrimaryQual".
 240                      " object to write_seq() as a parameter named \"source\"");
 241     }
 242     my $header = ($source->can("header") && $source->header) ?
 243         $source->header :
 244             ($source->can("id") && $source->id) ?
 245                 $source->id :
 246                     "unknown";
 247     my @quals = $source->qual();
 248     # ::dumpValue(\@quals);
 249     my $desc = $source->desc if $source->can('desc');
 250     $desc ||= '';
 251     $self->_print (">$header $desc\n");
 252     my (@slice,$max,$length);
 253     $length = $source->length();
 254
 255     if ( not(defined($oneline)) || $oneline == 0) {
 256         # $width quality values per line
 257         for (my $count = 1; $count<=$length; $count+= $width) {
 258             if ($count+$width > $length) {
 259                 $max = $length;
 260             } else {
 261                 $max = $count+$width-1;
 262             }
 263             my @slice = @{$source->subqual($count,$max)};
 264             $self->_print (join(' ',@slice), "\n");
 265         }
 266     } else {
 267         # quality values on a single line
 268         my @slice = @{$source->qual};
 269         $self->_print (join(' ',@slice), "\n");
 270     }
 271
 272     $self->flush if $self->_flush_on_write && defined $self->_fh;
 273     return 1;
 274 }
 275
 276
 277 1;
 278
 279 __END__