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 Maasha::Calc;
  33 use Data::Dumper;
  34
  35 use vars qw( @ISA @EXPORT_OK );
  36
  37 require Exporter;
  38
  39 @ISA = qw( Exporter );
  40
  41 use constant {
  42     SEQ_NAME => 0,
  43     SEQ      => 1,
  44     SCORE    => 2,
  45 };
  46
  47 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  48
  49
  50 sub solexa_get_entry_octal
  51 {
  52     # Martin A. Hansen, August 2008
  53
  54     # Get the next Solexa entry form a file and returns
  55     # a triple of [ seq_name, seq, score ].
  56     # We asume a Solexa entry consists of four lines:
  57
  58     # @HWI-EAS157_20FFGAAXX:2:1:888:434
  59     # TTGGTCGCTCGCTCCGCGACCTCAGATCAGACGTGG
  60     # +HWI-EAS157_20FFGAAXX:2:1:888:434
  61     # hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhKhe
  62
  63     my ( $fh,   # filehandle to Solexa file
  64        ) = @_;
  65
  66     # Returns a list.
  67
  68     my ( $seq_header, $seq, $score_head, $score );
  69
  70     $seq_header   = <$fh>;
  71     $seq          = <$fh>;
  72     $score_header = <$fh>;
  73     $score        = <$fh>;
  74
  75     return if not defined $score;
  76
  77     chomp $seq_header;
  78     chomp $seq;
  79     chomp $score_header;
  80     chomp $score;
  81
  82     $seq_header =~ s/^@//;
  83     $score      =~ s/(.)/int( score_oct2dec( $1 ) ) . ";"/ge;
  84
  85     return wantarray ? ( $seq_header, $seq, $score ) : [ $seq_header, $seq, $score ];
  86 }
  87
  88
  89 sub score_oct2dec
  90 {
  91     # Martin A. Hansen, August 2008
  92
  93     # Converts a Solexa score in octal format to decimal format.
  94
  95     # http://maq.sourceforge.net/fastq.shtml
  96
  97     my ( $char,   # octal score
  98        ) = @_;
  99
 100     # Returns a float
 101
 102     return 10 * log( 1 + 10 ** ( ( ord( $char ) - 64 ) / 10.0 ) ) / log( 10 );
 103 }
 104
 105
 106 sub solexa_get_entry_decimal
 107 {
 108     # Martin A. Hansen, August 2008
 109
 110     # Get the next Solexa entry form a file and returns
 111     # a triple of [ seq_name, seq, score ].
 112     # We asume a Solexa entry consists of four lines:
 113
 114     # @USI-EAS18_131_4_1_352_619
 115     # ATGGATGGGTTGGAGATGCCCTCTGTAGGCACCAT
 116     # +USI-EAS18_131_4_1_352_619
 117     # 40 40 40 40 40 40 40 40 40 40 40 25 25 40 40 34 40 40 40 40 32 39 40 36 34 19 40 19 30 16 21 11 21 18 11
 118
 119     my ( $fh,   # filehandle to Solexa file
 120        ) = @_;
 121
 122     # Returns a list.
 123
 124     my ( $seq_header, $seq, $score_head, $score );
 125
 126     $seq_header   = <$fh>;
 127     $seq          = <$fh>;
 128     $score_header = <$fh>;
 129     $score        = <$fh>;
 130
 131     return if not defined $score;
 132
 133     chomp $seq_header;
 134     chomp $seq;
 135     chomp $score_header;
 136     chomp $score;
 137
 138     $seq_header =~ s/^@//;
 139     $score      =~ s/ /;/g;
 140
 141     return wantarray ? ( $seq_header, $seq, $score ) : [ $seq_header, $seq, $score ];
 142 }
 143
 144
 145 sub solexa2biopiece
 146 {
 147     # Martin A. Hansen, September 2008.
 148
 149     # Converts a Solexa entry to a Biopiece record.
 150
 151     my ( $entry,     # Solexa entry
 152          $quality,   # Quality cutoff
 153        ) = @_;
 154
 155     # Returns a hashref.
 156
 157     my ( @seqs, @scores, $i, %record );
 158
 159     @seqs   = split //,  $entry->[ SEQ ];
 160     @scores = split /;/, $entry->[ SCORE ];
 161
 162     for ( $i = 0; $i < @scores; $i++ ) {
 163         $seqs[ $i ] = lc $seqs[ $i ] if $scores[ $i ] < $quality;
 164     }
 165
 166     $record{ "SEQ_NAME" }   = $entry->[ SEQ_NAME ];
 167     $record{ "SEQ" }        = $entry->[ SEQ ];
 168     $record{ "SCORES" }     = $entry->[ SCORE ];
 169     $record{ "SEQ_LEN" }    = scalar @seqs;
 170     $record{ "SCORE_MEAN" } = sprintf ( "%.2f", Maasha::Calc::mean( \@scores ) );
 171
 172     return wantarray ? %record : \%record;
 173 }
 174
 175
 176 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 177
 178
 179 1;