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 get_entry
 320 {
 321     # Martin A. Hansen, July 2009.
 322
 323     # Gets the next FASTQ entry from a given filehandle.
 324
 325     my ( $fh,   # filehandle
 326        ) = @_;
 327
 328     # Returns a list
 329
 330     my ( $seq, $seq_name, $qual, $qual_name );
 331
 332     $seq_name  = <$fh>;
 333     $seq       = <$fh>;
 334     $qual_name = <$fh>;
 335     $qual      = <$fh>;
 336
 337     return unless $seq;
 338
 339     chomp $seq;
 340     chomp $seq_name;
 341     chomp $qual;
 342     chomp $qual_name;
 343
 344     $seq_name =~ s/^@//;
 345
 346     return wantarray ? ( $seq_name, $seq, $qual ) : [ $seq_name, $seq, $qual ];
 347 }
 348
 349
 350 sub put_entry
 351 {
 352     # Martin A. Hansen, July 2009.
 353
 354     # Output a FASTQ entry to STDOUT or a filehandle.
 355
 356     my ( $entry,   # FASTQ entry
 357          $fh,      # filehandle - OPTIONAL
 358        ) = @_;
 359
 360     # Returns nothing.
 361
 362     $fh ||= \*STDOUT;
 363
 364     print $fh "@" . $entry->[ SEQ_NAME ] . "\n";
 365     print $fh $entry->[ SEQ ] . "\n";
 366     print $fh "+\n";
 367     print $fh $entry->[ SCORES ] . "\n";
 368 }
 369
 370
 371 sub fastq2biopiece
 372 {
 373     # Martin A. Hansen, July 2009.
 374
 375     # Converts a FASTQ entry to a Biopiece record, where
 376     # the FASTQ quality scores are converted to numerics.
 377
 378     my ( $entry,     # FASTQ entry,
 379        ) = @_;
 380
 381     # Returns a hash.
 382
 383     my ( $record );
 384
 385     $record->{ 'SEQ' }        = $entry->[ SEQ ];
 386     $record->{ 'SEQ_NAME' }   = $entry->[ SEQ_NAME ];
 387     $record->{ 'SCORES' }     = $entry->[ SCORES ];
 388     $record->{ 'SEQ_LEN' }    = length $entry->[ SEQ ];
 389
 390     return wantarray ? %{ $record } : $record;
 391 }
 392
 393
 394 sub biopiece2fastq
 395 {
 396     # Martin A. Hansen, July 2009.
 397
 398     # Converts a Biopiece record to a FASTQ entry.
 399
 400     my ( $record,   # Biopiece record
 401        ) = @_;
 402
 403     # Returns a list.
 404
 405     my ( $list );
 406
 407     if ( exists $record->{ 'SEQ' } and exists $record->{ 'SEQ_NAME' } and exists $record->{ 'SCORES' } )
 408     {
 409         $list->[ SEQ_NAME ] = $record->{ 'SEQ_NAME' };
 410         $list->[ SEQ ]      = $record->{ 'SEQ' };
 411         $list->[ SCORES ]   = $record->{ 'SCORES' };
 412
 413         $list->[ SCORES ] =~ s/(\d+);/chr( ( $1 <= 93 ? $1 : 93 ) + 33 )/ge;
 414
 415         return wantarray ? @{ $list } : $list;
 416     }
 417
 418     return;
 419 }
 420
 421
 422 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<