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, @scores );
  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     @scores = split //, $score;
  83
  84     map { $_ = int( score_oct2dec( $_ ) ) } @scores;
  85
  86     $seq_header =~ s/^@//;
  87
  88     $score = join( ";", @scores );
  89
  90     return wantarray ? ( $seq_header, $seq, $score ) : [ $seq_header, $seq, $score ];
  91 }
  92
  93
  94 sub score_oct2dec
  95 {
  96     # Martin A. Hansen, August 2008
  97
  98     # Converts a Solexa score in octal format to decimal format.
  99
 100     # http://maq.sourceforge.net/fastq.shtml
 101
 102     my ( $char,   # octal score
 103        ) = @_;
 104
 105     # Returns a float
 106
 107     return 10 * log(1 + 10 ** ((ord( $char ) - 64) / 10.0)) / log(10);
 108 }
 109
 110
 111 sub solexa_get_entry_decimal
 112 {
 113     # Martin A. Hansen, August 2008
 114
 115     # Get the next Solexa entry form a file and returns
 116     # a triple of [ seq_name, seq, score ].
 117     # We asume a Solexa entry consists of four lines:
 118
 119     # @USI-EAS18_131_4_1_352_619
 120     # ATGGATGGGTTGGAGATGCCCTCTGTAGGCACCAT
 121     # +USI-EAS18_131_4_1_352_619
 122     # 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
 123
 124     my ( $fh,   # filehandle to Solexa file
 125        ) = @_;
 126
 127     # Returns a list.
 128
 129     my ( $seq_header, $seq, $score_head, $score, @scores );
 130
 131     $seq_header   = <$fh>;
 132     $seq          = <$fh>;
 133     $score_header = <$fh>;
 134     $score        = <$fh>;
 135
 136     return if not defined $score;
 137
 138     chomp $seq_header;
 139     chomp $seq;
 140     chomp $score_header;
 141     chomp $score;
 142
 143     @scores = split / /, $score;
 144
 145     $seq_header =~ s/^@//;
 146
 147     $score = join( ";", @scores );
 148
 149     return wantarray ? ( $seq_header, $seq, $score ) : [ $seq_header, $seq, $score ];
 150 }
 151
 152
 153 sub solexa2biopiece
 154 {
 155     # Martin A. Hansen, September 2008.
 156
 157     # Converts a Solexa entry to a Biopiece record.
 158
 159     my ( $entry,     # Solexa entry
 160          $quality,   # Quality cutoff
 161        ) = @_;
 162
 163     # Returns a hashref.
 164
 165     my ( @seqs, @scores, $i, %record );
 166
 167     @seqs   = split //,  $entry->[ SEQ ];
 168     @scores = split /;/, $entry->[ SCORE ];
 169
 170     for ( $i = 0; $i < @scores; $i++ ) {
 171         $seqs[ $i ] = lc $seqs[ $i ] if $scores[ $i ] < $quality;
 172     }
 173
 174     $record{ "SEQ_NAME" }   = $entry->[ SEQ_NAME ];
 175     $record{ "SEQ" }        = $entry->[ SEQ ];
 176     $record{ "SCORES" }     = $entry->[ SCORE ];
 177     $record{ "SEQ_LEN" }    = scalar @seqs;
 178     $record{ "SCORE_MEAN" } = sprintf ( "%.2f", Maasha::Calc::mean( \@scores ) );
 179
 180     return wantarray ? %record : \%record;
 181 }
 182
 183
 184 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 185
 186
 187 1;