X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bp_bin%2Fplot_scores;h=915dd34aa3e42475cc8155ceab444b2b1b5f1d3a;hb=48bea5c28b89dc5586d0bddb338ccd6ba23aa1f9;hp=af37cc2a4639db430114b30ebd3429213c93a29b;hpb=c99b48ef600e53e67be86bd408a9efbd3204859b;p=biopieces.git

diff --git a/bp_bin/plot_scores b/bp_bin/plot_scores
index af37cc2..915dd34 100755
--- 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
@@ -20,75 +20,88 @@
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-# This program is part of the Biopieces framework (www.biopieces.org).
-
-use warnings;
-use strict;
-use Maasha::Biopieces;
-use Maasha::Fastq;
-use Maasha::Plot;
-
 # Plot a histogram of mean sequence quality scores.
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-
-require 'biopieces'
+require 'maasha/biopieces'
+require 'maasha/seq'
 require 'gnuplot'
-require 'pp'
+require 'narray'
 
 terminals = "dumb,x11,aqua,post,pdf,png,svg"
+title     = "Mean Quality Scores"
+xlabel    = "Sequence position"
+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}
+casts << {:long=>'xlabel',    :short=>'X', :type=>'string', :mandatory=>false, :default=>xlabel, :allowed=>nil,       :disallowed=>nil}
+casts << {:long=>'ylabel',    :short=>'Y', :type=>'string', :mandatory=>false, :default=>ylabel, :allowed=>nil,       :disallowed=>nil}
+
+options = Biopieces.options_parse(ARGV, casts)
+
+SCORES_MAX = 10_000
 
-bp = Biopieces.new
+scores_vec = NArray.int(SCORES_MAX)
+count_vec  = NArray.int(SCORES_MAX)
+max        = 0
 
-options = bp.parse(ARGV, casts)
+Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
+  input.each_record do |record|
+    if record[:SCORES]
+      scores = record[:SCORES]
 
-BASE_SOLEXA = 64
+      if scores.length > 0
+        raise BiopiecesError, "score string too long: #{scores.length} > #{SCORES_MAX}" if scores.length > SCORES_MAX
 
-sum_hash   = Hash.new(0)
-count_hash = Hash.new(0)
+        scores_vec[0 ... scores.length] += NArray.to_na(scores, "byte") - Seq::SCORE_BASE
+        count_vec[0 ... scores.length]  += 1
 
-bp.each_record do |record|
-  if record[:SCORES]
-    scores = record[:SCORES]
-    (0 ... scores.length).each do |i|
-      sum_hash[i]   += (scores[i].ord - BASE_SOLEXA)
-      count_hash[i] += 1
+        max = scores.length if scores.length > max
+      end
     end
-  end
 
-  bp.puts record unless options[:no_stream]
+    output.puts record unless options[:no_stream]
+  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 *= (Seq::SCORE_MAX / count_vec.max(0).to_f)
 
-(0 ... sum_hash.size).each do |i|
-  x << i
-  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|
     plot.terminal options[:terminal]
-    plot.title    "Mean Quality Scores"
-    plot.ylabel   "Mean score"
-    plot.xlabel   "Sequence position"
+    plot.title    options[:title]
+    plot.xlabel   options[:xlabel]
+    plot.ylabel   options[:ylabel]
     plot.output   options[:data_out] if options[:data_out]
     plot.xrange   "[#{x.min - 1}:#{x.max + 1}]"
-    plot.yrange   "[0:40]"
+    plot.yrange   "[#{Seq::SCORE_MIN}:#{Seq::SCORE_MAX}]"
     plot.style    "fill solid 0.5 border"
     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 lt rgb \"black\""
+        ds.title = "relative count"
+      end
     end
   end
 end