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 void 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    found = 0;
 257     int    i     = 0;
 258     int    pos   = -1;
 259     double sum   = 0;
 260     double mean  = 0.0;
 261
 262     if ( window_size > strlen( scores ) )
 263     {
 264         fprintf( stderr, "ERROR: window_size > scores string: %d > %d\n\n", window_size, strlen(scores) );
 265         exit( 1 );
 266     }
 267
 268     /* ---- fill up window ---- */
 269
 270     for ( i = 0; i < window_size; i++ ) {
 271         sum += solexa2dec( scores[ i ] );
 272     }
 273
 274     mean = sum / window_size;
 275
 276     if ( mean <= min ) {
 277         found = 1;
 278         pos   = 0;
 279     }
 280
 281     /* --- scan the rest of the scores ---- */
 282
 283     while ( ! found && i < strlen( scores ) )
 284     {
 285         sum += solexa2dec( scores[ i ] );
 286         sum -= solexa2dec( scores[ i - window_size ] );
 287
 288         mean = ( mean < sum / window_size ) ? mean : sum / window_size;
 289
 290         // printf( "char->%c   score->%d   sum->%f   mean->%f\n", scores[i], solexa2dec(scores[i]),sum, mean);
 291
 292         i++;
 293
 294         if ( mean <= min ) {
 295             found = 1;
 296             pos   = i - window_size;
 297         }
 298     }
 299
 300     Inline_Stack_Vars;
 301     Inline_Stack_Reset;
 302
 303     Inline_Stack_Push( sv_2mortal( newSViv( mean ) ) );
 304     Inline_Stack_Push( sv_2mortal( newSViv( pos ) ) );
 305
 306     Inline_Stack_Done;
 307 }
 308
 309
 310 void softmask_solexa_str( char *seq, char *scores, int threshold )
 311 {
 312     /* Martin A. Hansen, July 2009 */
 313
 314     /* Given a sequence string and a score string (in Solexa range) */
 315     /* lowercases all residues where the decimal score is below a given threshold. */
 316
 317     int i = 0;
 318
 319     for ( i = 0; i < strlen( seq ); i++ )
 320     {
 321         if ( solexa2dec( scores[ i ] ) < threshold ) {
 322             seq[ i ] = tolower( seq[ i ] );
 323         }
 324     }
 325 }
 326
 327
 328 void softmask_phred_str( char *seq, char *scores, int threshold )
 329 {
 330     /* Martin A. Hansen, July 2009 */
 331
 332     /* Given a sequence string and a score string (in Phred range) */
 333     /* lowercases all residues where the decimal score is below a given threshold. */
 334
 335     int i = 0;
 336
 337     for ( i = 0; i < strlen( seq ); i++ )
 338     {
 339         if ( phred2dec( scores[ i ] ) < threshold ) {
 340             seq[ i ] = tolower( seq[ i ] );
 341         }
 342     }
 343 }
 344
 345
 346 int trim_left( char *scores, int min )
 347 {
 348     /* Martin A. Hansen, June 2010 */
 349
 350     /* Starting from the left in a score string, */
 351     /* locate the position when the score is above */
 352     /* a given min.*/
 353
 354     int pos = 0;
 355
 356     while ( pos < strlen( scores ) && solexa2dec( scores[ pos ] ) <= min ) {
 357         pos++;
 358     }
 359
 360     return pos;
 361 }
 362
 363
 364 int trim_right( char *scores, int min )
 365 {
 366     /* Martin A. Hansen, June 2010 */
 367
 368     /* Starting from the right in a score string, */
 369     /* locate the position when the score is above */
 370     /* a given min.*/
 371
 372     int len = strlen( scores );
 373     int pos = len;
 374
 375     while ( pos > 0 && solexa2dec( scores[ pos ] ) <= min ) {
 376         pos--;
 377     }
 378
 379     if ( pos == 0 ) {
 380         pos = len;
 381     }
 382
 383     return pos;
 384 }
 385
 386 END_C
 387
 388
 389 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 390
 391
 392 sub solexa_str2dec_str
 393 {
 394     # Martin A. Hansen, September 2009
 395
 396     # Converts a string of Solexa octal scores to a ; separated
 397     # string of decimal scores.
 398
 399     my ( $scores,   # Solexa scores
 400        ) = @_;
 401
 402     # Returns a string.
 403
 404     $scores =~ s/(.)/ solexa2dec( $1 ) . ";"/eg;
 405
 406     return $scores;
 407 }
 408
 409
 410 sub phred_str2dec_str
 411 {
 412     # Martin A. Hansen, September 2009
 413
 414     # Converts a string of Phred octal scores to a ; separated
 415     # string of decimal scores.
 416
 417     my ( $scores,   # Phred scores
 418        ) = @_;
 419
 420     # Returns a string.
 421
 422     $scores =~ s/(.)/ phred2dec( $1 ) . ";"/eg;
 423
 424     return $scores;
 425 }
 426
 427
 428 sub dec_str2solexa_str
 429 {
 430     # Martin A. Hansen, May 2010.
 431
 432     # Converts a ; separated string of decimal scores to a
 433     # string of Solexa octal scores.
 434
 435     my ( $scores,   # Decimal score string
 436        ) = @_;
 437
 438     # Returns a string.
 439
 440     $scores =~ s/(-\d{1,2})/0/g;
 441     $scores =~ s/(\d{1,2});?/dec2solexa( $1 )/eg;
 442
 443     return $scores;
 444 }
 445
 446 sub dec_str2phred_str
 447 {
 448     # Martin A. Hansen, November 2013.
 449
 450     # Converts a ; separated string of decimal scores to a
 451     # string of Phred scores.
 452
 453     my ( $scores,   # Decimal score string
 454        ) = @_;
 455
 456     # Returns a string.
 457
 458     $scores =~ s/(-\d{1,2})/0/g;
 459     $scores =~ s/(\d{1,2});?/dec2phred( $1 )/eg;
 460
 461     return $scores;
 462 }
 463
 464 sub get_entry
 465 {
 466     # Martin A. Hansen, July 2009.
 467
 468     # Gets the next FASTQ entry from a given filehandle.
 469
 470     my ( $fh,   # filehandle
 471        ) = @_;
 472
 473     # Returns a list
 474
 475     my ( $seq, $seq_name, $qual, $qual_name );
 476
 477     $seq_name  = <$fh>;
 478     $seq       = <$fh>;
 479     $qual_name = <$fh>;
 480     $qual      = <$fh>;
 481
 482     return unless $seq;
 483
 484     chomp $seq;
 485     chomp $seq_name;
 486     chomp $qual;
 487     chomp $qual_name;
 488
 489     $seq_name =~ s/^@//;
 490
 491     return wantarray ? ( $seq_name, $seq, $qual ) : [ $seq_name, $seq, $qual ];
 492 }
 493
 494
 495 sub put_entry
 496 {
 497     # Martin A. Hansen, July 2009.
 498
 499     # Output a FASTQ entry to STDOUT or a filehandle.
 500
 501     my ( $entry,   # FASTQ entry
 502          $fh,      # filehandle - OPTIONAL
 503        ) = @_;
 504
 505     # Returns nothing.
 506
 507     $fh ||= \*STDOUT;
 508
 509     print $fh "@" . $entry->[ SEQ_NAME ] . "\n";
 510     print $fh $entry->[ SEQ ] . "\n";
 511     print $fh "+\n";
 512     print $fh $entry->[ SCORES ] . "\n";
 513 }
 514
 515
 516 sub fastq2biopiece
 517 {
 518     # Martin A. Hansen, July 2009.
 519
 520     # Converts a FASTQ entry to a Biopiece record, where
 521     # the FASTQ quality scores are converted to numerics.
 522
 523     my ( $entry,     # FASTQ entry,
 524        ) = @_;
 525
 526     # Returns a hash.
 527
 528     my ( $record );
 529
 530     $record->{ 'SEQ' }        = $entry->[ SEQ ];
 531     $record->{ 'SEQ_NAME' }   = $entry->[ SEQ_NAME ];
 532     $record->{ 'SCORES' }     = $entry->[ SCORES ];
 533     $record->{ 'SEQ_LEN' }    = length $entry->[ SEQ ];
 534
 535     return wantarray ? %{ $record } : $record;
 536 }
 537
 538
 539 sub biopiece2fastq
 540 {
 541     # Martin A. Hansen, July 2009.
 542
 543     # Converts a Biopiece record to a FASTQ entry.
 544
 545     my ( $record,   # Biopiece record
 546        ) = @_;
 547
 548     # Returns a list.
 549
 550     my ( $list );
 551
 552     if ( exists $record->{ 'SEQ' } and exists $record->{ 'SEQ_NAME' } and exists $record->{ 'SCORES' } )
 553     {
 554         $list->[ SEQ_NAME ] = $record->{ 'SEQ_NAME' };
 555         $list->[ SEQ ]      = $record->{ 'SEQ' };
 556         $list->[ SCORES ]   = $record->{ 'SCORES' };
 557
 558         $list->[ SCORES ] =~ s/(\d+);/chr( ( $1 <= 93 ? $1 : 93 ) + 33 )/ge;
 559
 560         return wantarray ? @{ $list } : $list;
 561     }
 562
 563     return;
 564 }
 565
 566
 567 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<