options = Biopieces.options_parse(ARGV, casts)
options = Biopieces.options_parse(ARGV, casts)
-ILLUMINA_BASE = 64
-ILLUMINA_MIN = 0
-ILLUMINA_MAX = 40
-SCORES_MAX = 10_000
scores_vec = NArray.int(SCORES_MAX)
count_vec = NArray.int(SCORES_MAX)
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
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
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
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
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.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"
plot.style "fill solid 0.5 border"
plot.xtics "out"
plot.ytics "out"
# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
require 'maasha/biopieces'
require 'maasha/biopieces'
# Expading class Hash with possibly evil monkey patch.
class Hash
# Expading class Hash with possibly evil monkey patch.
class Hash
def scores2dec!
if self.has_key? :SCORES
self[:SCORES].gsub! /./ do |score|
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! ";"
end
self[:SCORES].chomp! ";"
run "$bp -I $in -O $tmp"
assert_no_diff $tmp $out.1
clean
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
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")
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
# Method to calculate a consensus quality from a Consensus object.
def consensus_qual
# 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
"GTG" => "V", "GCG" => "A", "GAG" => "E", "GGG" => "G"
}
"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
# 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
include Digest
include Trim
def quality_trim_right(min_qual, min_len = 1)
check_trim_args(min_qual)
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)
def quality_trim_right!(min_qual, min_len = 1)
check_trim_args(min_qual)
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)
def quality_trim_left(min_qual, min_len = 1)
check_trim_args(min_qual)
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)
def quality_trim_left!(min_qual, min_len = 1)
check_trim_args(min_qual)
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)
def quality_trim(min_qual, min_len = 1)
check_trim_args(min_qual)
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
pos_left = pos_right if pos_left > pos_right
def quality_trim!(min_qual, min_len = 1)
check_trim_args(min_qual)
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
pos_left = pos_right if pos_left > pos_right
def check_trim_args(min_qual)
raise TrimError, "no sequence" if self.seq.nil?
raise TrimError, "no quality score" if self.qual.nil?
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}"