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, tmpdir, file_pattern, file_patscan, cpus)
  41     @options      = options
  42     @file_pattern = file_pattern
  43     @file_patscan = file_patscan
  44     @cpus         = cpus
  45     @files_fasta  = Dir.glob(File.join(tmpdir, "*.fna"))
  46
  47     pat = Pattern.new(@options)
  48     pat.write(@file_pattern)
  49   end
  50
  51 #  def run
  52 #    child_count = 0
  53 #
  54 #    @files_fasta.each do |file|
  55 #      if fork
  56 #        Process.wait if ( child_count += 1 ) >= @cpus
  57 #      else
  58 #        command = command_compile(file)
  59 #        system(command)
  60 #        raise PatScanError, "Command failed: #{command}" unless $?.success?
  61 #        exit
  62 #      end
  63 #    end
  64 #  end
  65
  66   def run
  67     child_count = 0
  68
  69     @files_fasta.each do |file|
  70       Thread.pass while child_count >= @cpus
  71       child_count += 1
  72
  73       Thread.new do
  74         command = command_compile(file)
  75         system(command)
  76         raise PatScanError, "Command failed: #{command}" unless $?.success?
  77         child_count -= 1
  78       end
  79     end
  80   end
  81
  82 #  def run
  83 #    child_count = 0
  84 #    ch_mutex = Mutex.new
  85 #    threads = []
  86 #
  87 #    @files_fasta.each do |file|
  88 #      Thread.pass while child_count >= @cpus
  89 #      ch_mutex.synchronize { child_count += 1 }
  90 #
  91 #      threads << Thread.new do
  92 #        command = command_compile(file)
  93 #        system(command)
  94 #        raise PatScanError, "Command failed: #{command}" unless $?.success?
  95 #        ch_mutex.synchronize { child_count -= 1 }
  96 #      end
  97 #    end
  98 #
  99 #    threads.each { |t| t.join }
 100 #  end
 101
 102   def command_compile(file)
 103     commands = []
 104     commands << "nice -n 19"
 105     commands << "scan_for_matches"
 106     commands << @file_pattern
 107     commands << "< #{file}"
 108     commands << "> #{file}.out"
 109     command = commands.join(" ")
 110   end
 111
 112   def parse_results
 113     matches = {}
 114
 115     Fasta.open(@file_patscan, mode='r') do |ios|
 116       ios.each do |entry|
 117         if entry.seq_name =~ /^(\d+):\[(\d+),(\d+)\]$/
 118           name  = $1.to_i
 119           start = $2.to_i - 1
 120           stop  = $3.to_i - 1
 121           matches[name] = [start, stop - start + 1] unless matches.has_key? name
 122         else
 123           raise "Failed to parse sequence name: #{entry.seq_name}"
 124         end
 125       end
 126     end
 127
 128     matches
 129   end
 130 end
 131
 132 # Error class for Pattern errors.
 133 class PatternError < StandardError; end;
 134
 135 class Pattern
 136   def initialize(options)
 137     @options  = options
 138     @patterns = []
 139     @patterns << pattern_internal
 140     @patterns += patterns_end if @options[:partial]
 141   end
 142
 143   def to_i
 144     new_patterns = []
 145
 146     while @patterns.size > 1
 147       new_patterns = @patterns[0 ... -2]
 148       new_patterns << "( #{@patterns[-2 .. -1].join(' | ')} )"
 149
 150       @patterns = new_patterns
 151     end
 152
 153     @patterns.first
 154   end
 155
 156   def write(file)
 157     File.open(file, mode='w') do |ios|
 158       ios.puts self.to_i
 159     end
 160   end
 161
 162   private
 163
 164   def pattern_internal
 165     pattern = @options[:adaptor]
 166     mis = mis_count(pattern)
 167     ins = ins_count(pattern)
 168     del = del_count(pattern)
 169
 170     "#{pattern}[#{mis},#{ins},#{del}]"
 171   end
 172
 173   def patterns_end
 174     patterns = []
 175     adaptor  = @options[:adaptor]
 176
 177     raise PatternError, "len > adaptor length: #{@options[:len]} > #{adaptor.length - 1}" if @options[:len] > adaptor.length - 1
 178
 179     (adaptor.length - 1).downto(@options[:len]) do |i|
 180       pattern = adaptor[0 ... i]
 181       mis = mis_count(pattern)
 182       ins = ins_count(pattern)
 183       del = del_count(pattern)
 184       patterns << "#{pattern}[#{mis},#{ins},#{del}] $"
 185     end
 186
 187     patterns
 188   end
 189
 190   def mis_count(pattern)
 191     (pattern.length * @options[:mismatches] * 0.01).round
 192   end
 193
 194   def ins_count(pattern)
 195     (pattern.length * @options[:insertions] * 0.01).round
 196   end
 197
 198   def del_count(pattern)
 199     (pattern.length * @options[:deletions] * 0.01).round
 200   end
 201 end
 202
 203 casts = []
 204 casts << {:long=>'adaptor',    :short=>'a', :type=>'string', :mandatory=>true,  :default=>nil, :allowed=>nil, :disallowed=>nil}
 205 casts << {:long=>'partial',    :short=>'p', :type=>'flag',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 206 casts << {:long=>'len',        :short=>'l', :type=>'uint',   :mandatory=>false, :default=>10,  :allowed=>nil, :disallowed=>'0'}
 207 casts << {:long=>'mismatches', :short=>'m', :type=>'uint',   :mandatory=>false, :default=>10,  :allowed=>nil, :disallowed=>nil}
 208 casts << {:long=>'insertions', :short=>'i', :type=>'uint',   :mandatory=>false, :default=>5,   :allowed=>nil, :disallowed=>nil}
 209 casts << {:long=>'deletions',  :short=>'d', :type=>'uint',   :mandatory=>false, :default=>5,   :allowed=>nil, :disallowed=>nil}
 210 casts << {:long=>'cpus',       :short=>'c', :type=>'uint',   :mandatory=>false, :default=>1,   :allowed=>nil, :disallowed=>'0'}
 211
 212 BASE_PER_FILE = 10_000_000
 213
 214 options = Biopieces.options_parse(ARGV, casts)
 215
 216 #tmpdir       = Biopieces.mktmpdir
 217 tmpdir       = "Tyt" # DEBUG TODO
 218 file_records = File.join(tmpdir, "data.stream")
 219 file_pattern = File.join(tmpdir, "pattern.txt")
 220 file_patscan = File.join(tmpdir, "patscan.out")
 221
 222 number_file = 0
 223 number_seq  = 0
 224 bases       = 0
 225
 226 Biopieces.open(options[:stream_in], file_records) do |input, output|
 227   file_fasta = File.join(tmpdir, "#{number_file}.fna")
 228   out_fa     = Fasta.open(file_fasta, mode='w')
 229
 230   input.each do |record|
 231     output.puts record
 232
 233     if record.has_key? :SEQ
 234       record[:SEQ_NAME] = number_seq
 235       out_fa.puts record
 236
 237       number_seq += 1;
 238       bases      += record[:SEQ].length
 239
 240       if bases > BASE_PER_FILE
 241         out_fa.close
 242         bases = 0
 243         number_file += 1
 244         file_fasta = File.join(tmpdir, "#{number_file}.fna")
 245         out_fa     = Fasta.open(file_fasta, mode='w')
 246       end
 247     end
 248   end
 249
 250   out_fa.close if out_fa.respond_to? :close
 251 end
 252
 253 patscan = PatScan.new(options, tmpdir, file_pattern, file_patscan, options[:cpus])
 254 patscan.run
 255 exit
 256 matches = patscan.parse_results
 257
 258 number_seq = 0
 259
 260 Biopieces.open(file_records, options[:stream_out]) do |input, output|
 261   input.each_record do |record|
 262     if record.has_key? :SEQ
 263       if matches.has_key? number_seq
 264         record[:ADAPTOR_POS] = matches[number_seq].first
 265         record[:ADAPTOR_LEN] = matches[number_seq].last
 266       end
 267
 268       number_seq += 1;
 269     end
 270
 271     output.puts record
 272   end
 273 end
 274
 275
 276 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 277
 278
 279 __END__