From dfe922540177aeaad25c2da595b9bac1a2205c27 Mon Sep 17 00:00:00 2001 From: martinahansen Date: Mon, 25 Feb 2013 19:17:08 +0000 Subject: [PATCH] fixed SCORE_BASE global var git-svn-id: http://biopieces.googlecode.com/svn/trunk@2101 74ccb610-7750-0410-82ae-013aeee3265d --- bp_bin/plot_scores | 11 ++++------- bp_bin/scores_to_dec | 6 ++---- bp_test/test/test_find_SNPs | 4 ++++ code_ruby/lib/maasha/align.rb | 4 ++-- code_ruby/lib/maasha/seq.rb | 9 +++++---- code_ruby/lib/maasha/seq/trim.rb | 20 ++++++++++---------- 6 files changed, 27 insertions(+), 27 deletions(-) diff --git a/bp_bin/plot_scores b/bp_bin/plot_scores index bf0023a..30ffeb3 100755 --- a/bp_bin/plot_scores +++ b/bp_bin/plot_scores @@ -45,10 +45,7 @@ casts << {:long=>'ylabel', :short=>'Y', :type=>'string', :mandatory=>false, : 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) @@ -62,7 +59,7 @@ Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output| 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 @@ -76,7 +73,7 @@ end 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 @@ -90,7 +87,7 @@ Gnuplot.open do |gp| 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" diff --git a/bp_bin/scores_to_dec b/bp_bin/scores_to_dec index 56381a5..2d91c2b 100755 --- a/bp_bin/scores_to_dec +++ b/bp_bin/scores_to_dec @@ -28,10 +28,8 @@ # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - require 'maasha/biopieces' - -ILLUMINA_BASE = 64 +require 'maasha/seq' # Expading class Hash with possibly evil monkey patch. class Hash @@ -40,7 +38,7 @@ class Hash def scores2dec! if self.has_key? :SCORES self[:SCORES].gsub! /./ do |score| - score = (score.ord - ILLUMINA_BASE).to_s + ";" + score = (score.ord - Seq::SCORE_BASE).to_s + ";" end self[:SCORES].chomp! ";" diff --git a/bp_test/test/test_find_SNPs b/bp_test/test/test_find_SNPs index 61c5709..760e1b0 100755 --- a/bp_test/test/test_find_SNPs +++ b/bp_test/test/test_find_SNPs @@ -5,3 +5,7 @@ source "$BP_DIR/bp_test/lib/test.sh" run "$bp -I $in -O $tmp" assert_no_diff $tmp $out.1 clean + +run "$bp -I $in -c -O $tmp" +assert_no_diff $tmp $out.2 +clean diff --git a/code_ruby/lib/maasha/align.rb b/code_ruby/lib/maasha/align.rb index 2b63d91..ffa41c4 100755 --- a/code_ruby/lib/maasha/align.rb +++ b/code_ruby/lib/maasha/align.rb @@ -262,7 +262,7 @@ class Align def na_add_entries @entries.each_with_index do |entry, i| @na_seq[true, i] = NArray.to_na(entry.seq.downcase.tr(TR_NUC, TR_HEX), "byte") - @na_qual[true, i] = NArray.to_na(entry.qual, "byte") - SCORE_BASE if @has_qual + @na_qual[true, i] = NArray.to_na(entry.qual, "byte") - Seq::SCORE_BASE if @has_qual end end @@ -375,7 +375,7 @@ class Align # Method to calculate a consensus quality from a Consensus object. def consensus_qual - (@na_qual.mean(1).round.to_type("byte") + SCORE_BASE).to_s + (@na_qual.mean(1).round.to_type("byte") + Seq::SCORE_BASE).to_s end end end diff --git a/code_ruby/lib/maasha/seq.rb b/code_ruby/lib/maasha/seq.rb index 8833acb..577d51a 100644 --- a/code_ruby/lib/maasha/seq.rb +++ b/code_ruby/lib/maasha/seq.rb @@ -67,15 +67,16 @@ TRANS_TAB11 = { "GTG" => "V", "GCG" => "A", "GAG" => "E", "GGG" => "G" } -# Quality scores bases -SCORE_BASE = 64 -SCORE_MIN = 0 -SCORE_MAX = 40 # Error class for all exceptions to do with Seq. class SeqError < StandardError; end class Seq + # Quality scores bases + SCORE_BASE = 64 + SCORE_MIN = 0 + SCORE_MAX = 40 + include Digest include Trim diff --git a/code_ruby/lib/maasha/seq/trim.rb b/code_ruby/lib/maasha/seq/trim.rb index a11159b..8efabf9 100644 --- a/code_ruby/lib/maasha/seq/trim.rb +++ b/code_ruby/lib/maasha/seq/trim.rb @@ -34,7 +34,7 @@ module Trim def quality_trim_right(min_qual, min_len = 1) check_trim_args(min_qual) - pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE) + pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE) self.subseq(0, pos) end @@ -44,7 +44,7 @@ module Trim def quality_trim_right!(min_qual, min_len = 1) check_trim_args(min_qual) - pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE) + pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE) self.subseq!(0, pos) end @@ -54,7 +54,7 @@ module Trim def quality_trim_left(min_qual, min_len = 1) check_trim_args(min_qual) - pos = trim_left_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE) + pos = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE) self.subseq(pos) end @@ -64,7 +64,7 @@ module Trim def quality_trim_left!(min_qual, min_len = 1) check_trim_args(min_qual) - pos = trim_left_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE) + pos = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE) self.subseq!(pos) end @@ -74,8 +74,8 @@ module Trim def quality_trim(min_qual, min_len = 1) check_trim_args(min_qual) - pos_right = trim_right_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE) - pos_left = trim_left_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE) + pos_right = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE) + pos_left = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE) pos_left = pos_right if pos_left > pos_right @@ -87,8 +87,8 @@ module Trim def quality_trim!(min_qual, min_len = 1) check_trim_args(min_qual) - pos_right = trim_right_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE) - pos_left = trim_left_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE) + pos_right = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE) + pos_left = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE) pos_left = pos_right if pos_left > pos_right @@ -102,8 +102,8 @@ module Trim def check_trim_args(min_qual) raise TrimError, "no sequence" if self.seq.nil? raise TrimError, "no quality score" if self.qual.nil? - unless (SCORE_MIN .. SCORE_MAX).include? min_qual - raise TrimError, "minimum quality value: #{min_qual} out of range #{SCORE_MIN} .. #{SCORE_MAX}" + unless (Seq::SCORE_MIN .. Seq::SCORE_MAX).include? min_qual + raise TrimError, "minimum quality value: #{min_qual} out of range #{Seq::SCORE_MIN} .. #{Seq::SCORE_MAX}" end end -- 2.39.2