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 solexa_str2dec_str
 261 {
 262     # Martin A. Hansen, September 2009
 263
 264     # Converts a string of Solexa octal scores to a ; separated
 265     # string of decimal scores.
 266
 267     my ( $scores,   # Solexa scores
 268        ) = @_;
 269
 270     # Returns a string.
 271
 272     $scores =~ s/(.)/ solexa2dec( $1 ) . ";"/eg;
 273
 274     return $scores;
 275 }
 276
 277
 278 sub phred_str2dec_str
 279 {
 280     # Martin A. Hansen, September 2009
 281
 282     # Converts a string of Phred octal scores to a ; separated
 283     # string of decimal scores.
 284
 285     my ( $scores,   # Phred scores
 286        ) = @_;
 287
 288     # Returns a string.
 289
 290     $scores =~ s/(.)/ phred2dec( $1 ) . ";"/eg;
 291
 292     return $scores;
 293 }
 294
 295
 296 sub get_entry
 297 {
 298     # Martin A. Hansen, July 2009.
 299
 300     # Gets the next FASTQ entry from a given filehandle.
 301
 302     my ( $fh,   # filehandle
 303        ) = @_;
 304
 305     # Returns a list
 306
 307     my ( $seq, $seq_name, $qual, $qual_name );
 308
 309     $seq_name  = <$fh>;
 310     $seq       = <$fh>;
 311     $qual_name = <$fh>;
 312     $qual      = <$fh>;
 313
 314     return unless $seq;
 315
 316     chomp $seq;
 317     chomp $seq_name;
 318     chomp $qual;
 319     chomp $qual_name;
 320
 321     $seq_name =~ s/^@//;
 322
 323     return wantarray ? ( $seq_name, $seq, $qual ) : [ $seq_name, $seq, $qual ];
 324 }
 325
 326
 327 sub put_entry
 328 {
 329     # Martin A. Hansen, July 2009.
 330
 331     # Output a FASTQ entry to STDOUT or a filehandle.
 332
 333     my ( $entry,   # FASTQ entry
 334          $fh,      # filehandle - OPTIONAL
 335        ) = @_;
 336
 337     # Returns nothing.
 338
 339     $fh ||= \*STDOUT;
 340
 341     print $fh "@" . $entry->[ SEQ_NAME ] . "\n";
 342     print $fh $entry->[ SEQ ] . "\n";
 343     print $fh "+\n";
 344     print $fh $entry->[ SCORES ] . "\n";
 345 }
 346
 347
 348 sub fastq2biopiece
 349 {
 350     # Martin A. Hansen, July 2009.
 351
 352     # Converts a FASTQ entry to a Biopiece record, where
 353     # the FASTQ quality scores are converted to numerics.
 354
 355     my ( $entry,     # FASTQ entry,
 356        ) = @_;
 357
 358     # Returns a hash.
 359
 360     my ( $record );
 361
 362     $record->{ 'SEQ' }        = $entry->[ SEQ ];
 363     $record->{ 'SEQ_NAME' }   = $entry->[ SEQ_NAME ];
 364     $record->{ 'SCORES' }     = $entry->[ SCORES ];
 365     $record->{ 'SEQ_LEN' }    = length $entry->[ SEQ ];
 366
 367     return wantarray ? %{ $record } : $record;
 368 }
 369
 370
 371 sub biopiece2fastq
 372 {
 373     # Martin A. Hansen, July 2009.
 374
 375     # Converts a Biopiece record to a FASTQ entry.
 376
 377     my ( $record,   # Biopiece record
 378        ) = @_;
 379
 380     # Returns a list.
 381
 382     my ( $list );
 383
 384     if ( exists $record->{ 'SEQ' } and exists $record->{ 'SEQ_NAME' } and exists $record->{ 'SCORES' } )
 385     {
 386         $list->[ SEQ_NAME ] = $record->{ 'SEQ_NAME' };
 387         $list->[ SEQ ]      = $record->{ 'SEQ' };
 388         $list->[ SCORES ]   = $record->{ 'SCORES' };
 389
 390         $list->[ SCORES ] =~ s/(\d+);/chr( ( $1 <= 93 ? $1 : 93 ) + 33 )/ge;
 391
 392         return wantarray ? @{ $list } : $list;
 393     }
 394
 395     return;
 396 }
 397
 398
 399 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<