code_perl/Maasha/Solexa.pm

   1 package Maasha::Solexa;
   2
   3
   4 # Copyright (C) 2007-2008 Martin A. Hansen.
   5
   6 # This program is free software; you can redistribute it and/or
   7 # modify it under the terms of the GNU General Public License
   8 # as published by the Free Software Foundation; either version 2
   9 # of the License, or (at your option) any later version.
  10
  11 # This program is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15
  16 # You should have received a copy of the GNU General Public License
  17 # along with this program; if not, write to the Free Software
  18 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19
  20 # http://www.gnu.org/copyleft/gpl.html
  21
  22
  23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  24
  25
  26 # Routines for manipulation Solid sequence files with di-base encoding.
  27
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31
  32 use vars qw( @ISA @EXPORT_OK );
  33
  34 require Exporter;
  35
  36 @ISA = qw( Exporter );
  37
  38
  39 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  40
  41
  42 sub solexa_get_entry
  43 {
  44     # Martin A. Hansen, August 2008
  45
  46     # Get the next Solexa entry form a file and returns
  47     # a triple of [ seq_name, seq, score ].
  48     # We asume a Solexa entry consists of four lines:
  49
  50     # @HWI-EAS157_20FFGAAXX:2:1:888:434
  51     # TTGGTCGCTCGCTCCGCGACCTCAGATCAGACGTGG
  52     # +HWI-EAS157_20FFGAAXX:2:1:888:434
  53     # hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhKhe
  54
  55     my ( $fh,   # filehandle to Solexa file
  56        ) = @_;
  57
  58     # Returns a list.
  59
  60     my ( $seq_header, $seq, $score_head, $score, @scores );
  61
  62     $seq_header   = <$fh>;
  63     $seq          = <$fh>;
  64     $score_header = <$fh>;
  65     $score        = <$fh>;
  66
  67     return if not defined $score;
  68
  69     chomp $seq_header;
  70     chomp $seq;
  71     chomp $score_header;
  72     chomp $score;
  73
  74     @scores = split //, $score;
  75
  76     map { $_ = score_oct2dec( $_ ) } @scores;
  77
  78     $seq_header =~ s/^@//;
  79
  80     $score = join( ":", @scores );
  81
  82     return wantarray ? ( $seq_header, $seq, $score ) : [ $seq_header, $seq, $score ];
  83 }
  84
  85
  86 sub score_oct2dec
  87 {
  88     # Martin A. Hansen, August 2008
  89
  90     # Converts a Solexa score in octal format to decimal format.
  91
  92     # http://maq.sourceforge.net/fastq.shtml
  93
  94     my ( $char,   # octal score
  95        ) = @_;
  96
  97     # Returns a float
  98
  99     return 10 * log(1 + 10 ** ((ord( $char ) - 64) / 10.0)) / log(10);
 100 }
 101
 102
 103 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 104
 105
 106 1;