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 #define BASE_SOLEXA 64
  75 #define BASE_PHRED 33
  76
  77
  78 int phred2dec( char c )
  79 {
  80     /* Martin A. Hansen, July 2009 */
  81
  82     /* Converts a Phred score in octal (a char) to a decimal score, */
  83     /* which is returned. */
  84
  85     int score = 0;
  86
  87     score = ( int ) c - BASE_PHRED;
  88
  89     return score;
  90 }
  91
  92
  93 int solexa2dec( char c )
  94 {
  95     /* Martin A. Hansen, July 2009 */
  96
  97     /* Converts a Solexa score in octal (a char) to a decimal score, */
  98     /* which is returned. */
  99
 100     int score = 0;
 101
 102     score = ( int ) c - BASE_SOLEXA;
 103
 104     return score;
 105 }
 106
 107
 108 char dec2phred( int score )
 109 {
 110     /* Martin A. Hansen, July 2009 */
 111
 112     /* Converts a decimal score to an octal score (a char) in the Phred range, */
 113     /* which is returned. */
 114
 115     char c = 0;
 116
 117     c = ( char ) score + BASE_PHRED;
 118
 119     return c;
 120 }
 121
 122
 123 char dec2solexa( int score )
 124 {
 125     /* Martin A. Hansen, September 2009 */
 126
 127     /* Converts a decimal score to an octal score (a char) in the Solexa range, */
 128     /* which is returned. */
 129
 130     char c = 0;
 131
 132     c = ( char ) score + BASE_SOLEXA;
 133
 134     return c;
 135 }
 136
 137
 138 char solexa2phred( char score )
 139 {
 140     /* Martin A. Hansen, September 2009 */
 141
 142     /* Converts an octal score in the Solexa range to an octal score */
 143     /* in the Pred range, which is returned. */
 144
 145     int  dec = 0;
 146     char c   = 0;
 147
 148     dec = solexa2dec( score );
 149     c   = dec2phred( c );
 150
 151     return c;
 152 }
 153
 154
 155 char phred2solexa( char score )
 156 {
 157     /* Martin A. Hansen, September 2009 */
 158
 159     /* Converts an octal score in the Solexa range to an octal score */
 160     /* in the Pred range, which is returned. */
 161
 162     int  dec = 0;
 163     char c   = 0;
 164
 165     dec = phred2dec( score );
 166     c   = dec2solexa( c );
 167
 168     return c;
 169 }
 170
 171
 172 void solexa_str2phred_str( char *scores )
 173 {
 174     /* Martin A. Hansen, September 2009 */
 175
 176     /* Converts a string of Solexa octal scores to Phred octal scores in-line. */
 177
 178     int i = 0;
 179
 180     for ( i = 0; i < strlen( scores ); i++ ) {
 181         scores[ i ] = solexa2phred( scores[ i ] );
 182     }
 183 }
 184
 185
 186 void phred_str2solexa_str( char *scores )
 187 {
 188     /* Martin A. Hansen, September 2009 */
 189
 190     /* Converts a string of Phred octal scores to Solexa octal scores in-line. */
 191
 192     int i = 0;
 193
 194     for ( i = 0; i < strlen( scores ); i++ ) {
 195         scores[ i ] = phred2solexa( scores[ i ] );
 196     }
 197 }
 198
 199
 200 double phred_str_mean( char *scores )
 201 {
 202     /* Martin A. Hansen, November 2009 */
 203
 204     /* Calculates the mean score as a float which is retuned. */
 205
 206     int    len  = 0;
 207     int    i    = 0;
 208     int    sum  = 0;
 209     double mean = 0.0;
 210
 211     len = strlen( scores );
 212
 213     for ( i = 0; i < len; i++ ) {
 214         sum += phred2dec( scores[ i ] );
 215     }
 216
 217     mean = ( double ) sum / ( double ) len;
 218
 219     return mean;
 220 }
 221
 222
 223 double solexa_str_mean( char *scores )
 224 {
 225     /* Martin A. Hansen, November 2009 */
 226
 227     /* Calculates the mean score as a float which is retuned. */
 228
 229     int    len  = 0;
 230     int    i    = 0;
 231     int    sum  = 0;
 232     double mean = 0.0;
 233
 234     len = strlen( scores );
 235
 236     for ( i = 0; i < len; i++ ) {
 237         sum += solexa2dec( scores[ i ] );
 238     }
 239
 240     mean = ( double ) sum / ( double ) len;
 241
 242     return mean;
 243 }
 244
 245
 246 double solexa_str_mean_window( char *scores, int window_size, double min )
 247 {
 248     /* Martin A. Hansen, June 2010. */
 249
 250     /* Scans a score string by running a sliding window across   */
 251     /* the string and for each position calculate the mean score */
 252     /* for the window. Terminates and returns mean score if this */
 253     /* is lower than a given minimum otherwise the smallest mean */
 254     /* score is returned. */
 255
 256     int    i    = 0;
 257     double sum  = 0;
 258     double mean = 0.0;
 259
 260     if ( window_size > strlen( scores ) )
 261     {
 262         fprintf( stderr, "ERROR: window_size > scores string: %d > %d\n\n", window_size, strlen(scores) );
 263         exit( 1 );
 264     }
 265
 266     /* ---- fill up window ---- */
 267
 268     for ( i = 0; i < window_size; i++ ) {
 269         sum += solexa2dec( scores[ i ] );
 270     }
 271
 272     mean = sum / window_size;
 273
 274     if ( mean < min ) {
 275         return mean;
 276     }
 277
 278     /* --- scan the rest of the scores ---- */
 279
 280     while ( i < strlen( scores ) )
 281     {
 282         sum += solexa2dec( scores[ i ] );
 283         sum -= solexa2dec( scores[ i - window_size ] );
 284
 285         mean = ( mean < sum / window_size ) ? mean : sum / window_size;
 286
 287         // printf( "char->%c   score->%d   sum->%f   mean->%f\n", scores[i], solexa2dec(scores[i]),sum, mean);
 288
 289         if ( mean < min ) {
 290             return mean;
 291         }
 292
 293         i++;
 294     }
 295
 296     return mean;
 297 }
 298
 299
 300 void softmask_solexa_str( char *seq, char *scores, int threshold )
 301 {
 302     /* Martin A. Hansen, July 2009 */
 303
 304     /* Given a sequence string and a score string (in Solexa range) */
 305     /* lowercases all residues where the decimal score is below a given threshold. */
 306
 307     int i = 0;
 308
 309     for ( i = 0; i < strlen( seq ); i++ )
 310     {
 311         if ( solexa2dec( scores[ i ] ) < threshold ) {
 312             seq[ i ] = tolower( seq[ i ] );
 313         }
 314     }
 315 }
 316
 317
 318 void softmask_phred_str( char *seq, char *scores, int threshold )
 319 {
 320     /* Martin A. Hansen, July 2009 */
 321
 322     /* Given a sequence string and a score string (in Phred range) */
 323     /* lowercases all residues where the decimal score is below a given threshold. */
 324
 325     int i = 0;
 326
 327     for ( i = 0; i < strlen( seq ); i++ )
 328     {
 329         if ( phred2dec( scores[ i ] ) < threshold ) {
 330             seq[ i ] = tolower( seq[ i ] );
 331         }
 332     }
 333 }
 334
 335
 336 int trim_left( char *scores, int min )
 337 {
 338     /* Martin A. Hansen, June 2010 */
 339
 340     /* Starting from the left in a score string, */
 341     /* locate the position when the score is above */
 342     /* a given min.*/
 343
 344     int pos = 0;
 345
 346     while ( pos < strlen( scores ) && solexa2dec( scores[ pos ] ) <= min ) {
 347         pos++;
 348     }
 349
 350     return pos;
 351 }
 352
 353
 354 int trim_right( char *scores, int min )
 355 {
 356     /* Martin A. Hansen, June 2010 */
 357
 358     /* Starting from the right in a score string, */
 359     /* locate the position when the score is above */
 360     /* a given min.*/
 361
 362     int len = strlen( scores );
 363     int pos = len;
 364
 365     while ( pos > 0 && solexa2dec( scores[ pos ] ) <= min ) {
 366         pos--;
 367     }
 368
 369     if ( pos == 0 ) {
 370         pos = len;
 371     }
 372
 373     return pos;
 374 }
 375
 376 END_C
 377
 378
 379 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 380
 381
 382 sub solexa_str2dec_str
 383 {
 384     # Martin A. Hansen, September 2009
 385
 386     # Converts a string of Solexa octal scores to a ; separated
 387     # string of decimal scores.
 388
 389     my ( $scores,   # Solexa scores
 390        ) = @_;
 391
 392     # Returns a string.
 393
 394     $scores =~ s/(.)/ solexa2dec( $1 ) . ";"/eg;
 395
 396     return $scores;
 397 }
 398
 399
 400 sub phred_str2dec_str
 401 {
 402     # Martin A. Hansen, September 2009
 403
 404     # Converts a string of Phred octal scores to a ; separated
 405     # string of decimal scores.
 406
 407     my ( $scores,   # Phred scores
 408        ) = @_;
 409
 410     # Returns a string.
 411
 412     $scores =~ s/(.)/ phred2dec( $1 ) . ";"/eg;
 413
 414     return $scores;
 415 }
 416
 417
 418 sub dec_str2solexa_str
 419 {
 420     # Martin A. Hansen, May 2010.
 421
 422     # Converts a ; separated string of decimal scores to a
 423     # string of Solexa octal scores.
 424
 425     my ( $scores,   # Decimal score string
 426        ) = @_;
 427
 428     # Returns a string.
 429
 430     $scores =~ s/(\d{1,2});?/dec2solexa( $1 )/eg;
 431
 432     return $scores;
 433 }
 434
 435
 436 sub get_entry
 437 {
 438     # Martin A. Hansen, July 2009.
 439
 440     # Gets the next FASTQ entry from a given filehandle.
 441
 442     my ( $fh,   # filehandle
 443        ) = @_;
 444
 445     # Returns a list
 446
 447     my ( $seq, $seq_name, $qual, $qual_name );
 448
 449     $seq_name  = <$fh>;
 450     $seq       = <$fh>;
 451     $qual_name = <$fh>;
 452     $qual      = <$fh>;
 453
 454     return unless $seq;
 455
 456     chomp $seq;
 457     chomp $seq_name;
 458     chomp $qual;
 459     chomp $qual_name;
 460
 461     $seq_name =~ s/^@//;
 462
 463     return wantarray ? ( $seq_name, $seq, $qual ) : [ $seq_name, $seq, $qual ];
 464 }
 465
 466
 467 sub put_entry
 468 {
 469     # Martin A. Hansen, July 2009.
 470
 471     # Output a FASTQ entry to STDOUT or a filehandle.
 472
 473     my ( $entry,   # FASTQ entry
 474          $fh,      # filehandle - OPTIONAL
 475        ) = @_;
 476
 477     # Returns nothing.
 478
 479     $fh ||= \*STDOUT;
 480
 481     print $fh "@" . $entry->[ SEQ_NAME ] . "\n";
 482     print $fh $entry->[ SEQ ] . "\n";
 483     print $fh "+\n";
 484     print $fh $entry->[ SCORES ] . "\n";
 485 }
 486
 487
 488 sub fastq2biopiece
 489 {
 490     # Martin A. Hansen, July 2009.
 491
 492     # Converts a FASTQ entry to a Biopiece record, where
 493     # the FASTQ quality scores are converted to numerics.
 494
 495     my ( $entry,     # FASTQ entry,
 496        ) = @_;
 497
 498     # Returns a hash.
 499
 500     my ( $record );
 501
 502     $record->{ 'SEQ' }        = $entry->[ SEQ ];
 503     $record->{ 'SEQ_NAME' }   = $entry->[ SEQ_NAME ];
 504     $record->{ 'SCORES' }     = $entry->[ SCORES ];
 505     $record->{ 'SEQ_LEN' }    = length $entry->[ SEQ ];
 506
 507     return wantarray ? %{ $record } : $record;
 508 }
 509
 510
 511 sub biopiece2fastq
 512 {
 513     # Martin A. Hansen, July 2009.
 514
 515     # Converts a Biopiece record to a FASTQ entry.
 516
 517     my ( $record,   # Biopiece record
 518        ) = @_;
 519
 520     # Returns a list.
 521
 522     my ( $list );
 523
 524     if ( exists $record->{ 'SEQ' } and exists $record->{ 'SEQ_NAME' } and exists $record->{ 'SCORES' } )
 525     {
 526         $list->[ SEQ_NAME ] = $record->{ 'SEQ_NAME' };
 527         $list->[ SEQ ]      = $record->{ 'SEQ' };
 528         $list->[ SCORES ]   = $record->{ 'SCORES' };
 529
 530         $list->[ SCORES ] =~ s/(\d+);/chr( ( $1 <= 93 ? $1 : 93 ) + 33 )/ge;
 531
 532         return wantarray ? @{ $list } : $list;
 533     }
 534
 535     return;
 536 }
 537
 538
 539 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<