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                 @tmpdir       = tmpdir
  43     @file_pattern = file_pattern
  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     files_result = Dir.glob(File.join(@tmpdir, "*.out"))
 114
 115     matches = {}
 116
 117     files_result.each do |file|
 118       Fasta.open(file, 'r') do |ios|
 119         ios.each do |entry|
 120           if entry.seq_name =~ /^(\d+):\[(\d+),(\d+)\]$/
 121             name  = $1.to_i
 122             start = $2.to_i - 1
 123             stop  = $3.to_i - 1
 124             matches[name] = [start, stop - start + 1] unless matches.has_key? name
 125           else
 126             raise "Failed to parse sequence name: #{entry.seq_name}"
 127           end
 128         end
 129       end
 130     end
 131
 132     matches
 133   end
 134 end
 135
 136 # Error class for Pattern errors.
 137 class PatternError < StandardError; end;
 138
 139 class Pattern
 140   def initialize(options)
 141     @options  = options
 142     @patterns = []
 143     @patterns << pattern_internal
 144     @patterns += patterns_end if @options[:partial]
 145   end
 146
 147   def to_i
 148     new_patterns = []
 149
 150     while @patterns.size > 1
 151       new_patterns = @patterns[0 ... -2]
 152       new_patterns << "( #{@patterns[-2 .. -1].join(' | ')} )"
 153
 154       @patterns = new_patterns
 155     end
 156
 157     @patterns.first
 158   end
 159
 160   def write(file)
 161     File.open(file, 'w') do |ios|
 162       ios.puts self.to_i
 163     end
 164   end
 165
 166   private
 167
 168   def pattern_internal
 169     pattern = @options[:adaptor]
 170     mis = mis_count(pattern)
 171     ins = ins_count(pattern)
 172     del = del_count(pattern)
 173
 174     "#{pattern}[#{mis},#{ins},#{del}]"
 175   end
 176
 177   def patterns_end
 178     patterns = []
 179     adaptor  = @options[:adaptor]
 180
 181     raise PatternError, "len > adaptor length: #{@options[:len]} > #{adaptor.length - 1}" if @options[:len] > adaptor.length - 1
 182
 183     (adaptor.length - 1).downto(@options[:len]) do |i|
 184       pattern = adaptor[0 ... i]
 185       mis = mis_count(pattern)
 186       ins = ins_count(pattern)
 187       del = del_count(pattern)
 188       patterns << "#{pattern}[#{mis},#{ins},#{del}] $"
 189     end
 190
 191     patterns
 192   end
 193
 194   def mis_count(pattern)
 195     (pattern.length * @options[:mismatches] * 0.01).round
 196   end
 197
 198   def ins_count(pattern)
 199     (pattern.length * @options[:insertions] * 0.01).round
 200   end
 201
 202   def del_count(pattern)
 203     (pattern.length * @options[:deletions] * 0.01).round
 204   end
 205 end
 206
 207 casts = []
 208 casts << {:long=>'adaptor',    :short=>'a', :type=>'string', :mandatory=>true,  :default=>nil, :allowed=>nil, :disallowed=>nil}
 209 casts << {:long=>'partial',    :short=>'p', :type=>'flag',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 210 casts << {:long=>'len',        :short=>'l', :type=>'uint',   :mandatory=>false, :default=>10,  :allowed=>nil, :disallowed=>'0'}
 211 casts << {:long=>'mismatches', :short=>'m', :type=>'uint',   :mandatory=>false, :default=>10,  :allowed=>nil, :disallowed=>nil}
 212 casts << {:long=>'insertions', :short=>'i', :type=>'uint',   :mandatory=>false, :default=>5,   :allowed=>nil, :disallowed=>nil}
 213 casts << {:long=>'deletions',  :short=>'d', :type=>'uint',   :mandatory=>false, :default=>5,   :allowed=>nil, :disallowed=>nil}
 214 casts << {:long=>'cpus',       :short=>'c', :type=>'uint',   :mandatory=>false, :default=>1,   :allowed=>nil, :disallowed=>'0'}
 215
 216 BASE_PER_FILE = 10_000_000
 217
 218 options = Biopieces.options_parse(ARGV, casts)
 219
 220 tmpdir       = Biopieces.mktmpdir
 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, '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, '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__