X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bp_bin%2Ffind_SNPs;h=283bc26a68eb88ffd828832d0949ad3eb1ccce34;hb=6c1d90a31caa1fe54bedc965ba26940688a31ab0;hp=6671144cc3a70066e0a221f66e4a08ff9e638a54;hpb=004fa0ccbda503d45e8af18758706c261c13337a;p=biopieces.git

diff --git a/bp_bin/find_SNPs b/bp_bin/find_SNPs
index 6671144..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,27 +32,40 @@
 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|
@@ -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
 
-
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<