added coverage to find_SNPs

git-svn-id: http://biopieces.googlecode.com/svn/trunk@2100 74ccb610-7750-0410-82ae-013aeee3265d

diff --git a/bp_bin/find_SNPs b/bp_bin/find_SNPs
index 2b40f32d408d55abcea7b436015ec9b5c2ec6e9c..283bc26a68eb88ffd828832d0949ad3eb1ccce34 100755 (executable)
--- 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
 
-
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<