bp_bin/find_adaptor

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2012 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 # Remove adaptors or parts thereof from sequences in the stream.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31 require 'pp'
  32 require 'maasha/biopieces'
  33 require 'maasha/seq'
  34 require 'maasha/seq/backtrack'
  35
  36 def percent2real(length, percent)
  37   (length * percent * 0.01).round
  38 end
  39
  40 include BackTrack
  41
  42 casts = []
  43 casts << {:long=>'forward',     :short=>'f', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
  44 casts << {:long=>'forward_rc',  :short=>'F', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
  45 casts << {:long=>'reverse',     :short=>'r', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
  46 casts << {:long=>'reverse_rc',  :short=>'R', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
  47 casts << {:long=>'len_forward', :short=>'l', :type=>'uint',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>'0'}
  48 casts << {:long=>'len_reverse', :short=>'L', :type=>'uint',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>'0'}
  49 casts << {:long=>'mismatches',  :short=>'m', :type=>'uint',   :mandatory=>false, :default=>10,  :allowed=>nil, :disallowed=>nil}
  50 casts << {:long=>'insertions',  :short=>'i', :type=>'uint',   :mandatory=>false, :default=>5,   :allowed=>nil, :disallowed=>nil}
  51 casts << {:long=>'deletions',   :short=>'d', :type=>'uint',   :mandatory=>false, :default=>5,   :allowed=>nil, :disallowed=>nil}
  52
  53 options = Biopieces.options_parse(ARGV, casts)
  54
  55 if options[:forward_rc]
  56   options[:forward] = Seq.new("test", options[:forward_rc], 'dna').revcomp.seq
  57 end
  58
  59 if options[:reverse_rc]
  60   options[:reverse] = Seq.new("test", options[:reverse_rc], 'dna').revcomp.seq
  61 end
  62
  63 raise ArgumentError, "no adaptor specified" unless options[:forward] or options[:reverse]
  64
  65 if options[:forward]
  66   options[:len_forward] = options[:forward].length unless options[:len_forward]
  67
  68   if options[:len_forward] > options[:forward].length
  69     raise ArgumentError, "len_forward > forward adaptor (#{options[:len_forward]} > #{options[:forward].length})"
  70   end
  71
  72   fmis = percent2real(options[:forward].length, options[:mismatches])
  73   fins = percent2real(options[:forward].length, options[:insertions])
  74   fdel = percent2real(options[:forward].length, options[:deletions])
  75 end
  76
  77 if options[:reverse]
  78   options[:len_reverse] = options[:reverse].length unless options[:len_reverse]
  79
  80   if options[:len_reverse] > options[:reverse].length
  81     raise ArgumentError, "len_reverse > reverse adaptor (#{options[:len_reverse]} > #{options[:reverse].length})"
  82   end
  83
  84   rmis = percent2real(options[:reverse].length, options[:mismatches])
  85   rins = percent2real(options[:reverse].length, options[:insertions])
  86   rdel = percent2real(options[:reverse].length, options[:deletions])
  87 end
  88
  89 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
  90   input.each do |record|
  91     if record[:SEQ] and record[:SEQ].length > 0
  92       entry = Seq.new_bp(record)
  93
  94       if options[:forward]
  95         if m = entry.patmatch(options[:forward], 0, fmis, fins, fdel)
  96           record[:ADAPTOR_POS_LEFT] = m.pos
  97           record[:ADAPTOR_LEN_LEFT] = m.length
  98           record[:ADAPTOR_PAT_LEFT] = m.match
  99         elsif options[:len_forward] < options[:forward].length
 100           len = options[:forward].length - 1
 101           pat = options[:forward]
 102
 103           while len >= options[:len_forward]
 104             fmis = percent2real(len, options[:mismatches])
 105             fins = percent2real(len, options[:insertions])
 106             fdel = percent2real(len, options[:deletions])
 107
 108             pat = pat[1 ... pat.length]
 109 $stderr.puts pat
 110 $stderr.puts entry.seq[0 ... len]
 111             if m = entry.patmatch(pat, [0, len], fmis, fins, fdel)
 112               record[:ADAPTOR_POS_LEFT] = m.pos
 113               record[:ADAPTOR_LEN_LEFT] = m.length
 114               record[:ADAPTOR_PAT_LEFT] = m.match
 115
 116               break
 117             end
 118
 119             len -= 1
 120           end
 121         end
 122       end
 123
 124       if options[:reverse]
 125         if m = entry.patmatch(options[:reverse], 0, rmis, rins, rdel)
 126           record[:ADAPTOR_POS_RIGHT] = m.pos
 127           record[:ADAPTOR_LEN_RIGHT] = m.length
 128           record[:ADAPTOR_PAT_RIGHT] = m.match
 129         elsif options[:len_reverse] < options[:reverse].length
 130           len = options[:reverse].length - 1
 131           pat = options[:reverse]
 132
 133           while len >= options[:len_reverse]
 134             rmis = percent2real(len, options[:mismatches])
 135             rins = percent2real(len, options[:insertions])
 136             rdel = percent2real(len, options[:deletions])
 137
 138             pat = pat[0 ... pat.length - 1]
 139
 140             if m = entry.patmatch(pat, entry.length - len, rmis, rins, rdel)
 141               record[:ADAPTOR_POS_RIGHT] = m.pos
 142               record[:ADAPTOR_LEN_RIGHT] = m.length
 143               record[:ADAPTOR_PAT_RIGHT] = m.match
 144
 145               break
 146             end
 147
 148             len -= 1
 149           end
 150         end
 151       end
 152     end
 153
 154     output.puts record
 155   end
 156 end
 157
 158 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 159
 160
 161 __END__