1 package Maasha::Solexa;
4 # Copyright (C) 2007-2008 Martin A. Hansen.
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.
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.
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.
20 # http://www.gnu.org/copyleft/gpl.html
23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
26 # Routines for manipulation Solid sequence files with di-base encoding.
29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
35 use vars qw( @ISA @EXPORT_OK );
39 @ISA = qw( Exporter );
47 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
50 sub solexa_get_entry_octal
52 # Martin A. Hansen, August 2008
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:
58 # @HWI-EAS157_20FFGAAXX:2:1:888:434
59 # TTGGTCGCTCGCTCCGCGACCTCAGATCAGACGTGG
60 # +HWI-EAS157_20FFGAAXX:2:1:888:434
61 # hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhKhe
63 my ( $fh, # filehandle to Solexa file
68 my ( $seq_header, $seq, $score_head, $score );
72 $score_header = <$fh>;
75 return if not defined $score;
82 $seq_header =~ s/^@//;
83 $score =~ s/(.)/int( score_oct2dec( $1 ) ) . ";"/ge;
85 return wantarray ? ( $seq_header, $seq, $score ) : [ $seq_header, $seq, $score ];
91 # Martin A. Hansen, August 2008
93 # Converts a Solexa score in octal format to decimal format.
95 # http://maq.sourceforge.net/fastq.shtml
97 my ( $char, # octal score
102 return 10 * log( 1 + 10 ** ( ( ord( $char ) - 64 ) / 10.0 ) ) / log( 10 );
106 sub solexa_get_entry_decimal
108 # Martin A. Hansen, August 2008
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:
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
119 my ( $fh, # filehandle to Solexa file
124 my ( $seq_header, $seq, $score_head, $score );
128 $score_header = <$fh>;
131 return if not defined $score;
138 $seq_header =~ s/^@//;
141 return wantarray ? ( $seq_header, $seq, $score ) : [ $seq_header, $seq, $score ];
147 # Martin A. Hansen, September 2008.
149 # Converts a Solexa entry to a Biopiece record.
151 my ( $entry, # Solexa entry
152 $quality, # Quality cutoff
157 my ( @seqs, @scores, $i, %record );
159 @seqs = split //, $entry->[ SEQ ];
160 @scores = split /;/, $entry->[ SCORE ];
162 for ( $i = 0; $i < @scores; $i++ ) {
163 $seqs[ $i ] = lc $seqs[ $i ] if $scores[ $i ] < $quality;
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 ) );
172 return wantarray ? %record : \%record;
176 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<