added quality threshold to read_fastq

git-svn-id: http://biopieces.googlecode.com/svn/trunk@603 74ccb610-7750-0410-82ae-013aeee3265d

diff --git a/bp_bin/read_fastq b/bp_bin/read_fastq
index ce1db6d69ac2d277d988474c347930f4c8726fbf..8aeedf9e1ca52b6fe9ced2334ddee3a12379b4e4 100755 (executable)
--- a/bp_bin/read_fastq
+++ b/bp_bin/read_fastq
@@ -41,7 +41,8 @@ my ( $options, $in, $out, $record, $data_in, $num, $entry );
 $options = Maasha::Biopieces::parse_options(
     [
         { long => 'data_in', short => 'i', type => 'files!', mandatory => 'no', default => undef, allowed => undef, disallowed => undef },
-        { long => 'num',     short => 'n', type => 'uint',   mandatory => 'no', default => undef, allowed => undef, disallowed => '0' },
+        { long => 'num',     short => 'n', type => 'uint',   mandatory => 'no', default => undef, allowed => undef, disallowed => 0     },
+        { long => 'quality', short => 'q', type => 'uint',   mandatory => 'no', default => 20,    allowed => undef, disallowed => undef },
     ]   
 );
 
@@ -60,7 +61,10 @@ if ( $options->{ 'data_in' } )
 
     while ( $entry = Maasha::Fastq::get_entry( $data_in ) ) 
     {
-        if ( $record = Maasha::Fastq::fastq2biopiece( $entry ) ) {
+        if ( $record = Maasha::Fastq::fastq2biopiece( $entry ) )
+        {
+            lowercase_low_scores( $record, $options->{ 'quality' } );
+
             Maasha::Biopieces::put_record( $record, $out );
         }
 
@@ -79,6 +83,37 @@ Maasha::Biopieces::close_stream( $out );
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
 
+sub lowercase_low_scores
+{
+    # Martin A. Hansen, July 2009.
+
+    # Lowercase residues in sequence that are below a given score threshold.
+
+    # This one would be lovely to have written in C.
+
+    my ( $record,    # Biopiece record
+         $quality,   # quality threshold
+       ) = @_;
+
+    # Returns nothing.
+
+    my ( @seqs, @scores, $i );
+
+    @seqs   = split //, $record->{ 'SEQ' };
+    @scores = split /;/, $record->{ 'SCORES' };
+
+    for ( $i = 0; $i < @scores; $i++ ) {
+        $seqs[ $i ] = lc $seqs[ $i ] if $scores[ $i ] < $quality;
+    }
+
+    $record->{ 'SEQ' }    = join "",  @seqs;
+    $record->{ 'SCORES' } = join ";", @scores;
+}
+
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
 BEGIN
 {
     Maasha::Biopieces::status_set();

diff --git a/code_perl/Maasha/Fastq.pm b/code_perl/Maasha/Fastq.pm
index 514916613521328b3e57cfcb1f1524179c39c713..be3c87c33562ae525b9fc8b8b802cde4e6d386e0 100644 (file)
--- a/code_perl/Maasha/Fastq.pm
+++ b/code_perl/Maasha/Fastq.pm
@@ -107,7 +107,7 @@ sub fastq2biopiece
     # Converts a FASTQ entry to a Biopiece record, where
     # the FASTQ quality scores are converted to numerics.
 
-    my ( $entry,   # FASTQ entry,
+    my ( $entry,     # FASTQ entry,
        ) = @_;
 
     # Returns a hash.