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
 447 sub get_entry
 448 {
 449     # Martin A. Hansen, July 2009.
 450
 451     # Gets the next FASTQ entry from a given filehandle.
 452
 453     my ( $fh,   # filehandle
 454        ) = @_;
 455
 456     # Returns a list
 457
 458     my ( $seq, $seq_name, $qual, $qual_name );
 459
 460     $seq_name  = <$fh>;
 461     $seq       = <$fh>;
 462     $qual_name = <$fh>;
 463     $qual      = <$fh>;
 464
 465     return unless $seq;
 466
 467     chomp $seq;
 468     chomp $seq_name;
 469     chomp $qual;
 470     chomp $qual_name;
 471
 472     $seq_name =~ s/^@//;
 473
 474     return wantarray ? ( $seq_name, $seq, $qual ) : [ $seq_name, $seq, $qual ];
 475 }
 476
 477
 478 sub put_entry
 479 {
 480     # Martin A. Hansen, July 2009.
 481
 482     # Output a FASTQ entry to STDOUT or a filehandle.
 483
 484     my ( $entry,   # FASTQ entry
 485          $fh,      # filehandle - OPTIONAL
 486        ) = @_;
 487
 488     # Returns nothing.
 489
 490     $fh ||= \*STDOUT;
 491
 492     print $fh "@" . $entry->[ SEQ_NAME ] . "\n";
 493     print $fh $entry->[ SEQ ] . "\n";
 494     print $fh "+\n";
 495     print $fh $entry->[ SCORES ] . "\n";
 496 }
 497
 498
 499 sub fastq2biopiece
 500 {
 501     # Martin A. Hansen, July 2009.
 502
 503     # Converts a FASTQ entry to a Biopiece record, where
 504     # the FASTQ quality scores are converted to numerics.
 505
 506     my ( $entry,     # FASTQ entry,
 507        ) = @_;
 508
 509     # Returns a hash.
 510
 511     my ( $record );
 512
 513     $record->{ 'SEQ' }        = $entry->[ SEQ ];
 514     $record->{ 'SEQ_NAME' }   = $entry->[ SEQ_NAME ];
 515     $record->{ 'SCORES' }     = $entry->[ SCORES ];
 516     $record->{ 'SEQ_LEN' }    = length $entry->[ SEQ ];
 517
 518     return wantarray ? %{ $record } : $record;
 519 }
 520
 521
 522 sub biopiece2fastq
 523 {
 524     # Martin A. Hansen, July 2009.
 525
 526     # Converts a Biopiece record to a FASTQ entry.
 527
 528     my ( $record,   # Biopiece record
 529        ) = @_;
 530
 531     # Returns a list.
 532
 533     my ( $list );
 534
 535     if ( exists $record->{ 'SEQ' } and exists $record->{ 'SEQ_NAME' } and exists $record->{ 'SCORES' } )
 536     {
 537         $list->[ SEQ_NAME ] = $record->{ 'SEQ_NAME' };
 538         $list->[ SEQ ]      = $record->{ 'SEQ' };
 539         $list->[ SCORES ]   = $record->{ 'SCORES' };
 540
 541         $list->[ SCORES ] =~ s/(\d+);/chr( ( $1 <= 93 ? $1 : 93 ) + 33 )/ge;
 542
 543         return wantarray ? @{ $list } : $list;
 544     }
 545
 546     return;
 547 }
 548
 549
 550 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<