fixed encoding bug in read_454

[biopieces.git] / code_perl / Maasha / Fastq.pm

diff --git a/code_perl/Maasha/Fastq.pm b/code_perl/Maasha/Fastq.pm
index f65cee20de186c65bbd8bff96de018f8954b969d..d70bf3e6cbbf2675e47feec6523976ae92c7bc9e 100644 (file)
--- a/code_perl/Maasha/Fastq.pm
+++ b/code_perl/Maasha/Fastq.pm
@@ -71,6 +71,9 @@ See also:
 
 http://maq.sourceforge.net/fastq.shtml */
 
+#define BASE_SOLEXA 64
+#define BASE_PHRED 33
+
 
 int phred2dec( char c )
 {
@@ -81,7 +84,7 @@ int phred2dec( char c )
 
     int score = 0;
 
-    score = ( int ) c - 33;
+    score = ( int ) c - BASE_PHRED;
 
     return score;
 }
@@ -96,33 +99,12 @@ int solexa2dec( char c )
 
     int score = 0;
 
-    score = ( int ) c - 64;
+    score = ( int ) c - BASE_SOLEXA;
 
     return score;
 }
 
 
-// int solexa2dec( char c )
-// {
-//     /* Martin A. Hansen, July 2009 */
-// 
-//     /* Converts a Solexa score in octal (a char) to a decimal score, */
-//     /* which is returned. */
-// 
-//     /* http://maq.sourceforge.net/fastq.shtml */
-//     /* $Q = 10 * log(1 + 10 ** (ord($sq) - 64) / 10.0)) / log(10); */
-// 
-//     int score = 0;
-//     int ascii = ( int ) c - 64;
-// 
-//     score = 10 * log( 1 + pow( 10, ascii / 10 ) ) / log( 10 );
-// 
-// //    printf( "char: %c   ascii: %d   score: %d\n", c, ascii, score );
-// 
-//     return score;
-// }
-
-
 char dec2phred( int score )
 {
     /* Martin A. Hansen, July 2009 */
@@ -132,7 +114,7 @@ char dec2phred( int score )
 
     char c = 0;
 
-    c = ( char ) score + 33;
+    c = ( char ) score + BASE_PHRED;
 
     return c;
 }
@@ -147,7 +129,7 @@ char dec2solexa( int score )
 
     char c = 0;
 
-    c = ( char ) score + 64;
+    c = ( char ) score + BASE_SOLEXA;
 
     return c;
 }
@@ -238,6 +220,93 @@ double phred_str_mean( char *scores )
 }
 
 
+double solexa_str_mean( char *scores )
+{
+    /* Martin A. Hansen, November 2009 */
+
+    /* Calculates the mean score as a float which is retuned. */
+
+    int    len  = 0;
+    int    i    = 0;
+    int    sum  = 0;
+    double mean = 0.0;
+
+    len = strlen( scores );
+
+    for ( i = 0; i < len; i++ ) {
+        sum += solexa2dec( scores[ i ] );
+    }
+
+    mean = ( double ) sum / ( double ) len;
+
+    return mean;
+}
+
+
+void solexa_str_mean_window( char *scores, int window_size, double min )
+{
+    /* Martin A. Hansen, June 2010. */
+
+    /* Scans a score string by running a sliding window across   */
+    /* the string and for each position calculate the mean score */
+    /* for the window. Terminates and returns mean score if this */
+    /* is lower than a given minimum otherwise the smallest mean */
+    /* score is returned. */
+
+    int    found = 0;
+    int    i     = 0;
+    int    pos   = -1;
+    double sum   = 0;
+    double mean  = 0.0;
+
+    if ( window_size > strlen( scores ) )
+    {
+        fprintf( stderr, "ERROR: window_size > scores string: %d > %d\n\n", window_size, strlen(scores) );
+        exit( 1 );
+    }
+
+    /* ---- fill up window ---- */
+    
+    for ( i = 0; i < window_size; i++ ) {
+        sum += solexa2dec( scores[ i ] );
+    }
+
+    mean = sum / window_size;
+
+    if ( mean <= min ) {
+        found = 1;
+        pos   = 0;
+    }
+
+    /* --- scan the rest of the scores ---- */
+
+    while ( ! found && i < strlen( scores ) )
+    {
+        sum += solexa2dec( scores[ i ] );
+        sum -= solexa2dec( scores[ i - window_size ] );
+
+        mean = ( mean < sum / window_size ) ? mean : sum / window_size;
+
+        // printf( "char->%c   score->%d   sum->%f   mean->%f\n", scores[i], solexa2dec(scores[i]),sum, mean);
+
+        i++;
+
+        if ( mean <= min ) {
+            found = 1;
+            pos   = i - window_size;
+        }
+    }
+
+    Inline_Stack_Vars;
+    Inline_Stack_Reset;
+
+    Inline_Stack_Push( sv_2mortal( newSViv( mean ) ) );
+    Inline_Stack_Push( sv_2mortal( newSViv( pos ) ) );
+
+    Inline_Stack_Done;
+}
+
+
 void softmask_solexa_str( char *seq, char *scores, int threshold )
 {
     /* Martin A. Hansen, July 2009 */
@@ -274,6 +343,46 @@ void softmask_phred_str( char *seq, char *scores, int threshold )
 }
 
 
+int trim_left( char *scores, int min )
+{
+    /* Martin A. Hansen, June 2010 */
+
+    /* Starting from the left in a score string, */
+    /* locate the position when the score is above */
+    /* a given min.*/
+
+    int pos = 0;
+
+    while ( pos < strlen( scores ) && solexa2dec( scores[ pos ] ) <= min ) {
+        pos++;
+    }
+
+    return pos;
+}
+
+
+int trim_right( char *scores, int min )
+{
+    /* Martin A. Hansen, June 2010 */
+
+    /* Starting from the right in a score string, */
+    /* locate the position when the score is above */
+    /* a given min.*/
+
+    int len = strlen( scores );
+    int pos = len;
+
+    while ( pos > 0 && solexa2dec( scores[ pos ] ) <= min ) {
+        pos--;
+    }
+
+    if ( pos == 0 ) {
+        pos = len;
+    }
+
+    return pos;
+}
+
 END_C
 
 
@@ -328,11 +437,29 @@ sub dec_str2solexa_str
 
     # Returns a string.
 
+    $scores =~ s/(-\d{1,2})/0/g;
     $scores =~ s/(\d{1,2});?/dec2solexa( $1 )/eg;
 
     return $scores;
 }
 
+sub dec_str2phred_str
+{
+    # Martin A. Hansen, November 2013.
+
+    # Converts a ; separated string of decimal scores to a
+    # string of Phred scores.
+
+    my ( $scores,   # Decimal score string
+       ) = @_;
+
+    # Returns a string.
+
+    $scores =~ s/(-\d{1,2})/0/g;
+    $scores =~ s/(\d{1,2});?/dec2phred( $1 )/eg;
+
+    return $scores;
+}
 
 sub get_entry
 {