bp_bin/find_adaptor

   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 # Remove adaptors or parts thereof from sequences in the stream.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31
  32 require 'pp'
  33 require 'maasha/biopieces'
  34 require 'maasha/fasta'
  35
  36 # Error class for PatScan errors.
  37 class PatScanError < StandardError; end;
  38
  39 class PatScan
  40   def initialize(options, file_fasta, file_pattern, file_patscan)
  41     @options      = options
  42     @file_fasta   = file_fasta
  43     @file_pattern = file_pattern
  44     @file_patscan = file_patscan
  45
  46     pat = Pattern.new(@options)
  47     pat.write(@file_pattern)
  48   end
  49
  50   def run
  51     commands = []
  52     commands << "nice -n 19"
  53     commands << "scan_for_matches"
  54     commands << @file_pattern
  55     commands << "< #{@file_fasta}"
  56     commands << "> #{@file_patscan}"
  57     command = commands.join(" ")
  58     system(command)
  59     raise PatScanError, "Command failed: #{command}" unless $?.success?
  60   end
  61
  62   def parse_results
  63     matches = {}
  64
  65     Fasta.open(@file_patscan, mode='r') do |ios|
  66       ios.each do |entry|
  67         if entry.seq_name =~ /^(\d+):\[(\d+),(\d+)\]$/
  68           name  = $1.to_i
  69           start = $2.to_i - 1
  70           stop  = $3.to_i - 1
  71           matches[name] = [start, stop - start + 1] unless matches.has_key? name
  72         else
  73           raise "Failed to parse sequence name: #{entry.seq_name}"
  74         end
  75       end
  76     end
  77
  78     matches
  79   end
  80 end
  81
  82 # Error class for Pattern errors.
  83 class PatternError < StandardError; end;
  84
  85 class Pattern
  86   def initialize(options)
  87     @options  = options
  88     @patterns = []
  89     @patterns << pattern_internal
  90     @patterns += patterns_end if @options[:partial]
  91   end
  92
  93   def to_i
  94     new_patterns = []
  95
  96     while @patterns.size > 1
  97       new_patterns = @patterns[0 ... -2]
  98       new_patterns << "( #{@patterns[-2 .. -1].join(' | ')} )"
  99
 100       @patterns = new_patterns
 101     end
 102
 103     @patterns.first
 104   end
 105
 106   def write(file)
 107     File.open(file, mode='w') do |ios|
 108       ios.puts self.to_i
 109     end
 110   end
 111
 112   private
 113
 114   def pattern_internal
 115     pattern = @options[:adaptor]
 116     mis = mis_count(pattern)
 117     ins = ins_count(pattern)
 118     del = del_count(pattern)
 119
 120     "#{pattern}[#{mis},#{ins},#{del}]"
 121   end
 122
 123   def patterns_end
 124     patterns = []
 125     adaptor  = @options[:adaptor]
 126
 127     raise PatternError, "len > adaptor length: #{@options[:len]} > #{adaptor.length - 1}" if @options[:len] > adaptor.length - 1
 128
 129     (adaptor.length - 1).downto(@options[:len]) do |i|
 130       pattern = adaptor[0 ... i]
 131       mis = mis_count(pattern)
 132       ins = ins_count(pattern)
 133       del = del_count(pattern)
 134       patterns << "#{pattern}[#{mis},#{ins},#{del}] $"
 135     end
 136
 137     patterns
 138   end
 139
 140   def mis_count(pattern)
 141     (pattern.length * @options[:mismatches] * 0.01).round
 142   end
 143
 144   def ins_count(pattern)
 145     (pattern.length * @options[:insertions] * 0.01).round
 146   end
 147
 148   def del_count(pattern)
 149     (pattern.length * @options[:deletions] * 0.01).round
 150   end
 151 end
 152
 153 casts = []
 154 casts << {:long=>'adaptor',    :short=>'a', :type=>'string', :mandatory=>true,  :default=>nil, :allowed=>nil, :disallowed=>nil}
 155 casts << {:long=>'partial',    :short=>'p', :type=>'flag',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 156 casts << {:long=>'len',        :short=>'l', :type=>'uint',   :mandatory=>false, :default=>10,  :allowed=>nil, :disallowed=>'0'}
 157 casts << {:long=>'mismatches', :short=>'m', :type=>'uint',   :mandatory=>false, :default=>10,  :allowed=>nil, :disallowed=>nil}
 158 casts << {:long=>'insertions', :short=>'i', :type=>'uint',   :mandatory=>false, :default=>5,   :allowed=>nil, :disallowed=>nil}
 159 casts << {:long=>'deletions',  :short=>'d', :type=>'uint',   :mandatory=>false, :default=>5,   :allowed=>nil, :disallowed=>nil}
 160
 161 options = Biopieces.options_parse(ARGV, casts)
 162
 163 tmpdir       = Biopieces.mktmpdir
 164 file_fasta   = File.join(tmpdir, "data.fna")
 165 file_records = File.join(tmpdir, "data.stream")
 166 file_pattern = File.join(tmpdir, "pattern.txt")
 167 file_patscan = File.join(tmpdir, "patscan.fna")
 168
 169 count = 0
 170
 171 Biopieces.open(options[:stream_in], file_records) do |input, output|
 172   Fasta.open(file_fasta, mode='w') do |out_fa|
 173     input.each do |record|
 174       output.puts record
 175
 176       if record.has_key? :SEQ
 177         record[:SEQ_NAME] = count
 178         out_fa.puts record
 179
 180         count += 1;
 181       end
 182     end
 183   end
 184 end
 185
 186 patscan = PatScan.new(options, file_fasta, file_pattern, file_patscan)
 187 patscan.run
 188 matches = patscan.parse_results
 189
 190 count = 0
 191
 192 Biopieces.open(file_records, options[:stream_out]) do |input, output|
 193   input.each_record do |record|
 194     if record.has_key? :SEQ
 195       if matches.has_key? count
 196         record[:ADAPTOR_POS] = matches[count].first
 197         record[:ADAPTOR_LEN] = matches[count].last
 198       end
 199
 200       count += 1;
 201     end
 202
 203     output.puts record
 204   end
 205 end
 206
 207
 208 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 209
 210
 211 __END__