From 44f810b56c353bfe302469a8313b0aada285b27e Mon Sep 17 00:00:00 2001 From: martinahansen Date: Thu, 21 Feb 2013 18:57:55 +0000 Subject: [PATCH] added coverage to find_SNPs git-svn-id: http://biopieces.googlecode.com/svn/trunk@2100 74ccb610-7750-0410-82ae-013aeee3265d --- bp_bin/find_SNPs | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/bp_bin/find_SNPs b/bp_bin/find_SNPs index 2b40f32..283bc26 100755 --- a/bp_bin/find_SNPs +++ b/bp_bin/find_SNPs @@ -1,6 +1,6 @@ #!/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 @@ -32,17 +32,21 @@ 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| @@ -53,6 +57,15 @@ Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output| 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| @@ -61,7 +74,8 @@ Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output| 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 @@ -76,6 +90,7 @@ Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output| else type = 'MISMATCH' end + new_record[:TYPE] = type end @@ -84,7 +99,6 @@ Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output| end end - # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< -- 2.39.2