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