code_perl/Maasha/Fastq.pm

   1 package Maasha::Fastq;
   2
   3 # Copyright (C) 2006-2009 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Routines for manipulation of FASTQ files and FASTQ entries.
  26
  27 # http://maq.sourceforge.net/fastq.shtml
  28
  29
  30 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  31
  32
  33 use warnings;
  34 use strict;
  35 use Data::Dumper;
  36 use Maasha::Calc;
  37 use vars qw( @ISA @EXPORT );
  38
  39 use constant {
  40     SEQ_NAME => 0,
  41     SEQ      => 1,
  42     SCORES   => 2,
  43 };
  44
  45 @ISA = qw( Exporter );
  46
  47 use Inline ( C => <<'END_C', DIRECTORY => $ENV{ "BP_TMP" } );
  48
  49 /*
  50 About Phred or FASTQ format:
  51
  52 This format is an adaption of the fasta format that contains quality scores.
  53 However, the fastq format is not completely compatible with the fastq files
  54 currently in existence, which is read by various applications (for example,
  55 BioPerl). Because a larger dynamic range of quality scores is used, the
  56 quality scores are encoded in ASCII as 64+score, instead of the standard
  57 32+score. This method is used to avoid running into non-printable characters.
  58
  59
  60 About Solexa or SCARF format:
  61
  62 SCARF (Solexa compact ASCII read format)—This easy-to-parse text based format,
  63 stores all the information for a single read in one line. Allowed values are
  64 --numeric and --symbolic. --numeric outputs the quality values as a
  65 space-separated string of numbers. --symbolic outputs the quality values as a
  66 compact string of ASCII characters. Subtract 64 from the ASCII code to get the
  67 corresponding quality values. Under the current numbering scheme, quality
  68 values can theoretically be as low as -5, which has an ASCII encoding of 59=';'.
  69
  70 See also:
  71
  72 http://maq.sourceforge.net/fastq.shtml */
  73
  74
  75 int phred2dec( char c )
  76 {
  77     /* Martin A. Hansen, July 2009 */
  78
  79     /* Converts a Phred score in octal (a char) to a decimal score, */
  80     /* which is returned. */
  81
  82     int score = 0;
  83
  84     score = ( int ) c - 33;
  85
  86     return score;
  87 }
  88
  89
  90 int solexa2dec( char c )
  91 {
  92     /* Martin A. Hansen, July 2009 */
  93
  94     /* Converts a Solexa score in octal (a char) to a decimal score, */
  95     /* which is returned. */
  96
  97     int score = 0;
  98
  99     score = ( int ) c - 64;
 100
 101     return score;
 102 }
 103
 104
 105 // int solexa2dec( char c )
 106 // {
 107 //     /* Martin A. Hansen, July 2009 */
 108 //
 109 //     /* Converts a Solexa score in octal (a char) to a decimal score, */
 110 //     /* which is returned. */
 111 //
 112 //     /* http://maq.sourceforge.net/fastq.shtml */
 113 //     /* $Q = 10 * log(1 + 10 ** (ord($sq) - 64) / 10.0)) / log(10); */
 114 //
 115 //     int score = 0;
 116 //     int ascii = ( int ) c - 64;
 117 //
 118 //     score = 10 * log( 1 + pow( 10, ascii / 10 ) ) / log( 10 );
 119 //
 120 // //    printf( "char: %c   ascii: %d   score: %d\n", c, ascii, score );
 121 //
 122 //     return score;
 123 // }
 124
 125
 126 char dec2phred( int score )
 127 {
 128     /* Martin A. Hansen, July 2009 */
 129
 130     /* Converts a decimal score to an octal score (a char) in the Phred range, */
 131     /* which is returned. */
 132
 133     char c = 0;
 134
 135     c = ( char ) score + 33;
 136
 137     return c;
 138 }
 139
 140
 141 char dec2solexa( int score )
 142 {
 143     /* Martin A. Hansen, September 2009 */
 144
 145     /* Converts a decimal score to an octal score (a char) in the Solexa range, */
 146     /* which is returned. */
 147
 148     char c = 0;
 149
 150     c = ( char ) score + 64;
 151
 152     return c;
 153 }
 154
 155
 156 char solexa2phred( char score )
 157 {
 158     /* Martin A. Hansen, September 2009 */
 159
 160     /* Converts an octal score in the Solexa range to an octal score */
 161     /* in the Pred range, which is returned. */
 162
 163     int  dec = 0;
 164     char c   = 0;
 165
 166     dec = solexa2dec( score );
 167     c   = dec2phred( c );
 168
 169     return c;
 170 }
 171
 172
 173 char phred2solexa( char score )
 174 {
 175     /* Martin A. Hansen, September 2009 */
 176
 177     /* Converts an octal score in the Solexa range to an octal score */
 178     /* in the Pred range, which is returned. */
 179
 180     int  dec = 0;
 181     char c   = 0;
 182
 183     dec = phred2dec( score );
 184     c   = dec2solexa( c );
 185
 186     return c;
 187 }
 188
 189
 190 void solexa_str2phred_str( char *scores )
 191 {
 192     /* Martin A. Hansen, September 2009 */
 193
 194     /* Converts a string of Solexa octal scores to Phred octal scores in-line. */
 195
 196     int i = 0;
 197
 198     for ( i + 0; i < strlen( scores ); i++ ) {
 199         scores[ i ] = solexa2phred( scores[ i ] );
 200     }
 201 }
 202
 203
 204 void phred_str2solexa_str( char *scores )
 205 {
 206     /* Martin A. Hansen, September 2009 */
 207
 208     /* Converts a string of Phred octal scores to Solexa octal scores in-line. */
 209
 210     int i = 0;
 211
 212     for ( i + 0; i < strlen( scores ); i++ ) {
 213         scores[ i ] = phred2solexa( scores[ i ] );
 214     }
 215 }
 216
 217
 218 void softmask_solexa_str( char *seq, char *scores, int threshold )
 219 {
 220     /* Martin A. Hansen, July 2009 */
 221
 222     /* Given a sequence string and a score string (in Solexa range) */
 223     /* lowercases all residues where the decimal score is below a given threshold. */
 224
 225     int i = 0;
 226
 227     for ( i = 0; i < strlen( seq ); i++ )
 228     {
 229         if ( solexa2dec( scores[ i ] ) < threshold ) {
 230             seq[ i ] = tolower( seq[ i ] );
 231         }
 232     }
 233 }
 234
 235
 236 void softmask_phred_str( char *seq, char *scores, int threshold )
 237 {
 238     /* Martin A. Hansen, July 2009 */
 239
 240     /* Given a sequence string and a score string (in Phred range) */
 241     /* lowercases all residues where the decimal score is below a given threshold. */
 242
 243     int i = 0;
 244
 245     for ( i = 0; i < strlen( seq ); i++ )
 246     {
 247         if ( phred2dec( scores[ i ] ) < threshold ) {
 248             seq[ i ] = tolower( seq[ i ] );
 249         }
 250     }
 251 }
 252
 253
 254 END_C
 255
 256
 257 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 258
 259
 260 sub get_entry
 261 {
 262     # Martin A. Hansen, July 2009.
 263
 264     # Gets the next FASTQ entry from a given filehandle.
 265
 266     my ( $fh,   # filehandle
 267        ) = @_;
 268
 269     # Returns a list
 270
 271     my ( $seq, $seq_name, $qual, $qual_name );
 272
 273     $seq_name  = <$fh>;
 274     $seq       = <$fh>;
 275     $qual_name = <$fh>;
 276     $qual      = <$fh>;
 277
 278     return unless $seq;
 279
 280     chomp $seq;
 281     chomp $seq_name;
 282     chomp $qual;
 283     chomp $qual_name;
 284
 285     $seq_name =~ s/^@//;
 286
 287     return wantarray ? ( $seq_name, $seq, $qual ) : [ $seq_name, $seq, $qual ];
 288 }
 289
 290
 291 sub put_entry
 292 {
 293     # Martin A. Hansen, July 2009.
 294
 295     # Output a FASTQ entry to STDOUT or a filehandle.
 296
 297     my ( $entry,   # FASTQ entry
 298          $fh,      # filehandle - OPTIONAL
 299        ) = @_;
 300
 301     # Returns nothing.
 302
 303     $fh ||= \*STDOUT;
 304
 305     print $fh "@" . $entry->[ SEQ_NAME ] . "\n";
 306     print $fh $entry->[ SEQ ] . "\n";
 307     print $fh "+\n";
 308     print $fh $entry->[ SCORES ] . "\n";
 309 }
 310
 311
 312 sub fastq2biopiece
 313 {
 314     # Martin A. Hansen, July 2009.
 315
 316     # Converts a FASTQ entry to a Biopiece record, where
 317     # the FASTQ quality scores are converted to numerics.
 318
 319     my ( $entry,     # FASTQ entry,
 320        ) = @_;
 321
 322     # Returns a hash.
 323
 324     my ( $record );
 325
 326     $record->{ 'SEQ' }        = $entry->[ SEQ ];
 327     $record->{ 'SEQ_NAME' }   = $entry->[ SEQ_NAME ];
 328     $record->{ 'SCORES' }     = $entry->[ SCORES ];
 329     $record->{ 'SEQ_LEN' }    = length $entry->[ SEQ ];
 330
 331     return wantarray ? %{ $record } : $record;
 332 }
 333
 334
 335 sub biopiece2fastq
 336 {
 337     # Martin A. Hansen, July 2009.
 338
 339     # Converts a Biopiece record to a FASTQ entry.
 340
 341     my ( $record,   # Biopiece record
 342        ) = @_;
 343
 344     # Returns a list.
 345
 346     my ( $list );
 347
 348     if ( exists $record->{ 'SEQ' } and exists $record->{ 'SEQ_NAME' } and exists $record->{ 'SCORES' } )
 349     {
 350         $list->[ SEQ_NAME ] = $record->{ 'SEQ_NAME' };
 351         $list->[ SEQ ]      = $record->{ 'SEQ' };
 352         $list->[ SCORES ]   = $record->{ 'SCORES' };
 353
 354         $list->[ SCORES ] =~ s/(\d+);/chr( ( $1 <= 93 ? $1 : 93 ) + 33 )/ge;
 355
 356         return wantarray ? @{ $list } : $list;
 357     }
 358
 359     return;
 360 }
 361
 362
 363 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<