lib/Bio/SeqIO/qual.pm

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