code_perl/Maasha/Fastq.pm

   1 package Maasha::Fastq;
   2
   3 # Copyright (C) 2006-2009 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Routines for manipulation of FASTQ files and FASTQ entries.
  26
  27 # http://maq.sourceforge.net/fastq.shtml
  28
  29
  30 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  31
  32
  33 use warnings;
  34 use strict;
  35 use Data::Dumper;
  36 use Maasha::Calc;
  37 use vars qw( @ISA @EXPORT );
  38
  39 @ISA = qw( Exporter );
  40
  41 use constant {
  42     SEQ_NAME => 0,
  43     SEQ      => 1,
  44     SCORES   => 2,
  45 };
  46
  47
  48 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  49
  50
  51 sub get_entry
  52 {
  53     # Martin A. Hansen, July 2009.
  54
  55     # Gets the next FASTQ entry from a given filehandle.
  56
  57     my ( $fh,   # filehandle
  58        ) = @_;
  59
  60     # Returns a list
  61
  62     my ( $seq, $seq_name, $qual, $qual_name );
  63
  64     $seq_name  = <$fh>;
  65     $seq       = <$fh>;
  66     $qual_name = <$fh>;
  67     $qual      = <$fh>;
  68
  69     return unless $seq;
  70
  71     chomp $seq;
  72     chomp $seq_name;
  73     chomp $qual;
  74     chomp $qual_name;
  75
  76     $seq_name =~ s/^@//;
  77
  78     return wantarray ? ( $seq_name, $seq, $qual ) : [ $seq_name, $seq, $qual ];
  79 }
  80
  81
  82 sub put_entry
  83 {
  84     # Martin A. Hansen, July 2009.
  85
  86     # Output a FASTQ entry to STDOUT or a filehandle.
  87
  88     my ( $entry,   # FASTQ entry
  89          $fh,      # filehandle - OPTIONAL
  90        ) = @_;
  91
  92     # Returns nothing.
  93
  94     $fh ||= \*STDOUT;
  95
  96     print $fh "@" . $entry->[ SEQ_NAME ] . "\n";
  97     print $fh $entry->[ SEQ ] . "\n";
  98     print $fh "+\n";
  99     print $fh $entry->[ SCORES ] . "\n";
 100 }
 101
 102
 103 sub fastq2biopiece
 104 {
 105     # Martin A. Hansen, July 2009.
 106
 107     # Converts a FASTQ entry to a Biopiece record, where
 108     # the FASTQ quality scores are converted to numerics.
 109
 110     my ( $entry,     # FASTQ entry,
 111        ) = @_;
 112
 113     # Returns a hash.
 114
 115     my ( $record );
 116
 117     $record->{ 'SEQ' }        = $entry->[ SEQ ];
 118     $record->{ 'SEQ_NAME' }   = $entry->[ SEQ_NAME ];
 119     $record->{ 'SCORES' }     = $entry->[ SCORES ];
 120
 121     $record->{ 'SCORES' }     =~ s/(.)/ord( $1 ) - 33 . ";"/ge; # http://maq.sourceforge.net/fastq.shtml
 122     $record->{ 'SCORE_MEAN' } = sprintf( "%.2f", Maasha::Calc::mean( [ split /;/, $record->{ 'SCORES' } ] ) );
 123
 124     return wantarray ? %{ $record } : $record;
 125 }
 126
 127
 128 sub biopiece2fastq
 129 {
 130     # Martin A. Hansen, July 2009.
 131
 132     # Converts a Biopiece record to a FASTQ entry.
 133
 134     my ( $record,   # Biopiece record
 135        ) = @_;
 136
 137     # Returns a list.
 138
 139     my ( $list );
 140
 141     if ( exists $record->{ 'SEQ' } and exists $record->{ 'SEQ_NAME' } and exists $record->{ 'SCORES' } )
 142     {
 143         $list->[ SEQ_NAME ] = $record->{ 'SEQ_NAME' };
 144         $list->[ SEQ ]      = $record->{ 'SEQ' };
 145         $list->[ SCORES ]   = $record->{ 'SCORES' };
 146
 147         $list->[ SCORES ] =~ s/(\d+);/chr( ( $1 <= 93 ? $1 : 93 ) + 33 )/ge;
 148
 149         return wantarray ? @{ $list } : $list;
 150     }
 151
 152     return;
 153 }
 154
 155
 156 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<