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 double phred_str_mean( char *scores )
 219 {
 220     /* Martin A. Hansen, November 2009 */
 221
 222     /* Calculates the mean score as a float which is retuned. */
 223
 224     int    len  = 0;
 225     int    i    = 0;
 226     int    sum  = 0;
 227     double mean = 0.0;
 228
 229     len = strlen( scores );
 230
 231     for ( i = 0; i < len; i++ ) {
 232         sum += phred2dec( scores[ i ] );
 233     }
 234
 235     mean = ( double ) sum / ( double ) len;
 236
 237     return mean;
 238 }
 239
 240
 241 double solexa_str_mean( char *scores )
 242 {
 243     /* Martin A. Hansen, November 2009 */
 244
 245     /* Calculates the mean score as a float which is retuned. */
 246
 247     int    len  = 0;
 248     int    i    = 0;
 249     int    sum  = 0;
 250     double mean = 0.0;
 251
 252     len = strlen( scores );
 253
 254     for ( i = 0; i < len; i++ ) {
 255         sum += solexa2dec( scores[ i ] );
 256     }
 257
 258     mean = ( double ) sum / ( double ) len;
 259
 260     return mean;
 261 }
 262
 263
 264 void softmask_solexa_str( char *seq, char *scores, int threshold )
 265 {
 266     /* Martin A. Hansen, July 2009 */
 267
 268     /* Given a sequence string and a score string (in Solexa range) */
 269     /* lowercases all residues where the decimal score is below a given threshold. */
 270
 271     int i = 0;
 272
 273     for ( i = 0; i < strlen( seq ); i++ )
 274     {
 275         if ( solexa2dec( scores[ i ] ) < threshold ) {
 276             seq[ i ] = tolower( seq[ i ] );
 277         }
 278     }
 279 }
 280
 281
 282 void softmask_phred_str( char *seq, char *scores, int threshold )
 283 {
 284     /* Martin A. Hansen, July 2009 */
 285
 286     /* Given a sequence string and a score string (in Phred range) */
 287     /* lowercases all residues where the decimal score is below a given threshold. */
 288
 289     int i = 0;
 290
 291     for ( i = 0; i < strlen( seq ); i++ )
 292     {
 293         if ( phred2dec( scores[ i ] ) < threshold ) {
 294             seq[ i ] = tolower( seq[ i ] );
 295         }
 296     }
 297 }
 298
 299
 300 END_C
 301
 302
 303 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 304
 305
 306 sub solexa_str2dec_str
 307 {
 308     # Martin A. Hansen, September 2009
 309
 310     # Converts a string of Solexa octal scores to a ; separated
 311     # string of decimal scores.
 312
 313     my ( $scores,   # Solexa scores
 314        ) = @_;
 315
 316     # Returns a string.
 317
 318     $scores =~ s/(.)/ solexa2dec( $1 ) . ";"/eg;
 319
 320     return $scores;
 321 }
 322
 323
 324 sub phred_str2dec_str
 325 {
 326     # Martin A. Hansen, September 2009
 327
 328     # Converts a string of Phred octal scores to a ; separated
 329     # string of decimal scores.
 330
 331     my ( $scores,   # Phred scores
 332        ) = @_;
 333
 334     # Returns a string.
 335
 336     $scores =~ s/(.)/ phred2dec( $1 ) . ";"/eg;
 337
 338     return $scores;
 339 }
 340
 341
 342 sub dec_str2solexa_str
 343 {
 344     # Martin A. Hansen, May 2010.
 345
 346     # Converts a ; separated string of decimal scores to a
 347     # string of Solexa octal scores.
 348
 349     my ( $scores,   # Decimal score string
 350        ) = @_;
 351
 352     # Returns a string.
 353
 354     $scores =~ s/(\d{1,2});?/dec2solexa( $1 )/eg;
 355
 356     return $scores;
 357 }
 358
 359
 360 sub get_entry
 361 {
 362     # Martin A. Hansen, July 2009.
 363
 364     # Gets the next FASTQ entry from a given filehandle.
 365
 366     my ( $fh,   # filehandle
 367        ) = @_;
 368
 369     # Returns a list
 370
 371     my ( $seq, $seq_name, $qual, $qual_name );
 372
 373     $seq_name  = <$fh>;
 374     $seq       = <$fh>;
 375     $qual_name = <$fh>;
 376     $qual      = <$fh>;
 377
 378     return unless $seq;
 379
 380     chomp $seq;
 381     chomp $seq_name;
 382     chomp $qual;
 383     chomp $qual_name;
 384
 385     $seq_name =~ s/^@//;
 386
 387     return wantarray ? ( $seq_name, $seq, $qual ) : [ $seq_name, $seq, $qual ];
 388 }
 389
 390
 391 sub put_entry
 392 {
 393     # Martin A. Hansen, July 2009.
 394
 395     # Output a FASTQ entry to STDOUT or a filehandle.
 396
 397     my ( $entry,   # FASTQ entry
 398          $fh,      # filehandle - OPTIONAL
 399        ) = @_;
 400
 401     # Returns nothing.
 402
 403     $fh ||= \*STDOUT;
 404
 405     print $fh "@" . $entry->[ SEQ_NAME ] . "\n";
 406     print $fh $entry->[ SEQ ] . "\n";
 407     print $fh "+\n";
 408     print $fh $entry->[ SCORES ] . "\n";
 409 }
 410
 411
 412 sub fastq2biopiece
 413 {
 414     # Martin A. Hansen, July 2009.
 415
 416     # Converts a FASTQ entry to a Biopiece record, where
 417     # the FASTQ quality scores are converted to numerics.
 418
 419     my ( $entry,     # FASTQ entry,
 420        ) = @_;
 421
 422     # Returns a hash.
 423
 424     my ( $record );
 425
 426     $record->{ 'SEQ' }        = $entry->[ SEQ ];
 427     $record->{ 'SEQ_NAME' }   = $entry->[ SEQ_NAME ];
 428     $record->{ 'SCORES' }     = $entry->[ SCORES ];
 429     $record->{ 'SEQ_LEN' }    = length $entry->[ SEQ ];
 430
 431     return wantarray ? %{ $record } : $record;
 432 }
 433
 434
 435 sub biopiece2fastq
 436 {
 437     # Martin A. Hansen, July 2009.
 438
 439     # Converts a Biopiece record to a FASTQ entry.
 440
 441     my ( $record,   # Biopiece record
 442        ) = @_;
 443
 444     # Returns a list.
 445
 446     my ( $list );
 447
 448     if ( exists $record->{ 'SEQ' } and exists $record->{ 'SEQ_NAME' } and exists $record->{ 'SCORES' } )
 449     {
 450         $list->[ SEQ_NAME ] = $record->{ 'SEQ_NAME' };
 451         $list->[ SEQ ]      = $record->{ 'SEQ' };
 452         $list->[ SCORES ]   = $record->{ 'SCORES' };
 453
 454         $list->[ SCORES ] =~ s/(\d+);/chr( ( $1 <= 93 ? $1 : 93 ) + 33 )/ge;
 455
 456         return wantarray ? @{ $list } : $list;
 457     }
 458
 459     return;
 460 }
 461
 462
 463 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<