added relative count to plot_scores + rewrite for speed

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

diff --git a/bp_bin/plot_scores b/bp_bin/plot_scores
index 20d12e3dffab519884f4642224b78d0dade8ef70..ff476f47c7359ac4b1593063bc84e508929c2a09 100755 (executable)
--- a/bp_bin/plot_scores
+++ b/bp_bin/plot_scores
@@ -1,6 +1,6 @@
 #!/usr/bin/env ruby
 
-# Copyright (C) 2007-2011 Martin A. Hansen.
+# Copyright (C) 2007-2012 Martin A. Hansen.
 
 # This program is free software; you can redistribute it and/or
 # modify it under the terms of the GNU General Public License
@@ -24,9 +24,9 @@
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-
 require 'maasha/biopieces'
 require 'gnuplot'
+require 'narray'
 require 'pp'
 
 terminals = "dumb,x11,aqua,post,pdf,png,svg"
@@ -36,6 +36,7 @@ ylabel    = "Mean score"
 
 casts = []
 casts << {:long=>'no_stream', :short=>'x', :type=>'flag',   :mandatory=>false, :default=>nil,    :allowed=>nil,       :disallowed=>nil}
+casts << {:long=>'count',     :short=>'c', :type=>'flag',   :mandatory=>false, :default=>nil,    :allowed=>nil,       :disallowed=>nil}
 casts << {:long=>'data_out',  :short=>'o', :type=>'file',   :mandatory=>false, :default=>nil,    :allowed=>nil,       :disallowed=>nil}
 casts << {:long=>'terminal',  :short=>'t', :type=>'string', :mandatory=>false, :default=>'dumb', :allowed=>terminals, :disallowed=>nil}
 casts << {:long=>'title',     :short=>'T', :type=>'string', :mandatory=>false, :default=>title,  :allowed=>nil,       :disallowed=>nil}
@@ -47,17 +48,24 @@ options = Biopieces.options_parse(ARGV, casts)
 ILLUMINA_BASE = 64
 ILLUMINA_MIN  = 0
 ILLUMINA_MAX  = 40
+SCORES_MAX    = 10_000
 
-sum_hash   = Hash.new(0)
-count_hash = Hash.new(0)
+scores_vec = NArray.int(SCORES_MAX)
+count_vec  = NArray.int(SCORES_MAX)
+max        = 0
 
 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
   input.each_record do |record|
     if record[:SCORES]
       scores = record[:SCORES]
-      (0 ... scores.length).each do |i|
-        sum_hash[i]   += (scores[i].ord - ILLUMINA_BASE)
-        count_hash[i] += 1
+
+      if scores.length > 0
+        raise BiopiecesError, "score string too long: #{scores.length} > #{SCORES_MAX}" if scores.length > SCORES_MAX
+
+        scores_vec[0 ... scores.length] += NArray.to_na(scores, "byte") - ILLUMINA_BASE
+        count_vec[0 ... scores.length]  += 1
+
+        max = scores.length if scores.length > max
       end
     end
 
@@ -65,13 +73,14 @@ Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
   end
 end
 
-x = []
-y = []
+mean_vec   = NArray.sfloat(max)
+mean_vec   = scores_vec[0 ... max].to_f / count_vec[0 ... max]
+count_vec  = count_vec[0 ... max].to_f
+count_vec *= (ILLUMINA_MAX / count_vec.max(0).to_f)
 
-(0 ... sum_hash.size).each do |i|
-  x << i + 1
-  y << sum_hash[i].to_f / count_hash[i].to_f
-end
+x  = (1 .. max).to_a
+y1 = mean_vec.to_a
+y2 = count_vec.to_a
 
 Gnuplot.open do |gp|
   Gnuplot::Plot.new(gp) do |plot|
@@ -86,9 +95,16 @@ Gnuplot.open do |gp|
     plot.xtics    "out"
     plot.ytics    "out"
     
-    plot.data << Gnuplot::DataSet.new([x, y]) do |ds|
-      ds.with = "boxes"
-      ds.notitle
+    plot.data << Gnuplot::DataSet.new([x, y1]) do |ds|
+      ds.with  = "boxes"
+      ds.title = "mean score"
+    end
+
+    if options[:count]
+      plot.data << Gnuplot::DataSet.new([x, y2]) do |ds|
+        ds.with  = "lines"
+        ds.title = "relative count"
+      end
     end
   end
 end