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