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