#!/usr/bin/env ruby
-# Copyright (C) 2007-2011 Martin A. Hansen.
+# Copyright (C) 2007-2013 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
require 'maasha/biopieces'
require 'pp'
-options = Biopieces.options_parse(ARGV)
+casts = []
+casts << {:long=>'coverage', :short=>'c', :type=>'flag', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
-snp_hash = Hash.new { |h,k| h[k] = Hash.new { |h,k| h[k] = Hash.new(0) } }
+options = Biopieces.options_parse(ARGV, casts)
+
+snp_hash = Hash.new { |hash, key| hash[key] = Hash.new { |h, k| h[k] = Hash.new(0) } }
+cov_hash = Hash.new { |hash, key| hash[key] = Hash.new(0) }
Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
input.each_record do |record|
output.puts record
- if record.has_key? :ALIGN and
- record.has_key? :S_ID and
- record.has_key? :S_BEG and
+ if record[:ALIGN] and
+ record[:S_ID] and
+ record[:S_BEG] and
record[:ALIGN] != '.'
record[:ALIGN].split(',').each do |snp|
pos, event = snp.split(':')
- pos += record[:S_BEG]
+ pos = pos.to_i + record[:S_BEG].to_i
- snp_hash[record[:S_ID].to_sym][pos.to_sym][event.to_sym] += 1
+ snp_hash[record[:S_ID].to_sym][pos][event.to_sym] += 1
end
end
+
+ if options[:coverage] and
+ record[:S_ID] and
+ record[:S_BEG] and
+ record[:S_END]
+ (record[:S_BEG].to_i .. record[:S_END].to_i).each do |i|
+ cov_hash[record[:S_ID].to_sym][i] += 1
+ end
+ end
end
snp_hash.each_pair do |s_id, pos_hash|
new_record[:S_ID] = s_id
pos_hash.each_pair do |pos, event_hash|
- new_record[:POS] = pos
+ new_record[:POS] = pos
+ new_record[:COVERAGE] = cov_hash[s_id][pos] if options[:coverage]
event_hash.each_pair do |event, count|
new_record[:EVENT] = event
else
type = 'MISMATCH'
end
+
new_record[:TYPE] = type
end
end
end
-
# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<