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