bp_bin/find_barcodes

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2011 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  22
  23 # This program is part of the Biopieces framework (www.biopieces.org).
  24
  25 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  26
  27 # Find barcodes in sequences in the stream.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31 require 'maasha/biopieces'
  32 require 'maasha/roche'
  33 require 'maasha/bits'
  34 require 'pp'
  35
  36 class BarCodeFinderError < StandardError; end
  37
  38 class BarCodeFinder
  39   def initialize(pos, max_mismatches)
  40     @pos            = pos
  41     @max_mismatches = max_mismatches
  42     @barcode_hash   = {}
  43     @barcode_sizes  = []
  44   end
  45
  46   def parse_barcodes(file)
  47     length = 0
  48     hash   = {}
  49
  50     File.open(file, "r") do |ios|
  51       while not ios.eof?
  52         name, barcode = ios.readline.chomp.split(/\s+/)
  53
  54         if length > 0
  55           raise BarCodeFinderError, "Uneven barcode length detected" if length != barcode.length
  56         else
  57           length = barcode.length
  58         end
  59
  60         hash[barcode.upcase.to_sym] = name
  61       end
  62     end
  63
  64     load_barcodes(hash)
  65   end
  66
  67   def load_barcodes(hash)
  68     @barcode_sizes << hash.keys.first.to_s.size
  69
  70     @barcode_hash.merge!(hash)
  71   end
  72
  73   def find_barcode(seq)
  74     raise BarCodeFinderError, "No barcodes to find" if @barcode_hash.empty?
  75
  76     @barcode_sizes.each do |size|
  77       hamming_dist = 0
  78       barcode      = seq[@pos ... @pos + size].upcase.to_sym
  79
  80       if @barcode_hash[barcode]
  81         return BarCode.new(barcode, @barcode_hash[barcode], @pos, size, hamming_dist)
  82       elsif @max_mismatches > 0
  83         @barcode_hash.each_key do |key|
  84           if key.to_s.length == barcode.to_s.size
  85             hamming_dist = barcode.to_s.hamming_distance(key.to_s)
  86
  87             if hamming_dist <= @max_mismatches
  88               return BarCode.new(key, @barcode_hash[key], @pos, size, hamming_dist)
  89             end
  90           end
  91         end
  92       end
  93     end
  94
  95     nil
  96   end
  97
  98   private
  99
 100   # Class to hold a BacCode object.
 101   class BarCode
 102     attr_accessor :barcode, :name, :pos, :len, :mismatches
 103
 104     def initialize(barcode, name, pos, len, mismatches)
 105       @barcode    = barcode
 106       @name       = name
 107       @pos        = pos
 108       @len        = len
 109       @mismatches = mismatches
 110     end
 111
 112     def to_hash
 113       hash = {}
 114       hash[:BARCODE]            = @barcode
 115       hash[:BARCODE_NAME]       = @name
 116       hash[:BARCODE_POS]        = @pos
 117       hash[:BARCODE_LEN]        = @len
 118       hash[:BARCODE_MISMATCHES] = @mismatches
 119       hash
 120     end
 121
 122     def start
 123       @pos + @len
 124     end
 125   end
 126 end
 127
 128 casts = []
 129 casts << {:long=>'barcodes_in', :short=>'b', :type=>'file!', :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 130 casts << {:long=>'pos',         :short=>'p', :type=>'uint',  :mandatory=>false, :default=>0,   :allowed=>nil,     :disallowed=>nil}
 131 casts << {:long=>'mismatches',  :short=>'m', :type=>'uint',  :mandatory=>false, :default=>0,   :allowed=>"0,1,2", :disallowed=>nil}
 132 casts << {:long=>'gsmids',      :short=>'g', :type=>'flag',  :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 133 casts << {:long=>'rlmids',      :short=>'r', :type=>'flag',  :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 134 casts << {:long=>'remove',      :short=>'R', :type=>'flag',  :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 135
 136 options = Biopieces.options_parse(ARGV, casts)
 137
 138 bc_finder = BarCodeFinder.new(options[:pos], options[:mismatches])
 139 bc_finder.parse_barcodes(options[:barcodes_in]) if options[:barcodes_in]
 140 bc_finder.load_barcodes(Roche::GSMID_HASH)      if options[:gsmids]
 141 bc_finder.load_barcodes(Roche::RLMID_HASH)      if options[:rlmids]
 142
 143 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
 144   input.each_record do |record|
 145     if record[:SEQ]
 146       if barcode = bc_finder.find_barcode(record[:SEQ])
 147         record.merge!(barcode.to_hash)
 148
 149         if options[:remove]
 150           record[:SEQ]     = record[:SEQ][barcode.start .. -1]
 151           record[:SCORES]  = record[:SCORES][barcode.start .. -1] if record[:SCORES]
 152           record[:SEQ_LEN] = record[:SEQ].length
 153         end
 154       end
 155     end
 156
 157     output.puts record
 158   end
 159 end
 160
 161
 162 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 163
 164
 165 __END__