# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
require 'maasha/biopieces'
+require 'maasha/seq'
require 'gnuplot'
require 'narray'
-require 'pp'
terminals = "dumb,x11,aqua,post,pdf,png,svg"
title = "Mean Quality Scores"
options = Biopieces.options_parse(ARGV, casts)
-ILLUMINA_BASE = 64
-ILLUMINA_MIN = 0
-ILLUMINA_MAX = 40
-SCORES_MAX = 10_000
+SCORES_MAX = 10_000
scores_vec = NArray.int(SCORES_MAX)
count_vec = NArray.int(SCORES_MAX)
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
+ scores_vec[0 ... scores.length] += NArray.to_na(scores, "byte") - Seq::SCORE_BASE
count_vec[0 ... scores.length] += 1
max = scores.length if scores.length > max
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)
+count_vec *= (Seq::SCORE_MAX / count_vec.max(0).to_f)
x = (1 .. max).to_a
y1 = mean_vec.to_a
plot.ylabel options[:ylabel]
plot.output options[:data_out] if options[:data_out]
plot.xrange "[#{x.min - 1}:#{x.max + 1}]"
- plot.yrange "[#{ILLUMINA_MIN}:#{ILLUMINA_MAX}]"
+ plot.yrange "[#{Seq::SCORE_MIN}:#{Seq::SCORE_MAX}]"
plot.style "fill solid 0.5 border"
plot.xtics "out"
plot.ytics "out"
if options[:count]
plot.data << Gnuplot::DataSet.new([x, y2]) do |ds|
- ds.with = "lines"
+ ds.with = "lines lt rgb \"black\""
ds.title = "relative count"
end
end