bp_bin/find_adaptor

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2012 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(tmpdir, cpus)
  41                 @tmpdir       = tmpdir
  42     @cpus         = cpus
  43     @files_fasta  = Dir.glob(File.join(@tmpdir, "*.fna"))
  44   end
  45
  46   def run(pattern)
  47     file_pat = File.join(@tmpdir, "pat.txt")
  48
  49     File.open(file_pat, 'w') do |ios|
  50       ios.puts pattern
  51     end
  52
  53     child_count = 0
  54     ch_mutex    = Mutex.new
  55     threads     = []
  56
  57     @files_fasta.each do |file_fasta|
  58       Thread.pass while child_count >= @cpus
  59       ch_mutex.synchronize { child_count += 1 }
  60
  61       threads << Thread.new do
  62         command = command_compile(file_pat, file_fasta)
  63         system(command)
  64         raise PatScanError, "Command failed: #{command}" unless $?.success?
  65         ch_mutex.synchronize { child_count -= 1 }
  66       end
  67     end
  68
  69     threads.each { |t| t.join }
  70   end
  71
  72   def parse_results
  73     files_result = Dir.glob(File.join(@tmpdir, "*.out"))
  74
  75     matches = {}
  76
  77     files_result.each do |file|
  78       Fasta.open(file, 'r') do |ios|
  79         ios.each do |entry|
  80           if entry.seq_name =~ /^(\d+):\[(\d+),(\d+)\]$/
  81             name  = $1.to_i
  82             start = $2.to_i - 1
  83             stop  = $3.to_i - 1
  84             matches[name] = [start, stop - start + 1] unless matches[name]
  85           else
  86             raise "Failed to parse sequence name: #{entry.seq_name}"
  87           end
  88         end
  89       end
  90     end
  91
  92     matches
  93   end
  94
  95   def clean
  96     Dir.glob(File.join(@tmpdir, "*.out")).each do |file|
  97       File.unlink(file)
  98     end
  99   end
 100
 101   private
 102
 103   def command_compile(file_pat, file_fasta)
 104     commands = []
 105     commands << "nice -n 19"
 106     commands << "scan_for_matches"
 107     commands << file_pat
 108     commands << "< #{file_fasta}"
 109     commands << "> #{file_fasta}.out"
 110     command = commands.join(" ")
 111   end
 112 end
 113
 114 # Error class for Pattern errors.
 115 class PatternError < StandardError; end;
 116
 117 class Pattern
 118   def self.create_left(adaptor, misp, insp, delp, len)
 119     patterns = []
 120
 121     (adaptor.length).downto(len) do |i|
 122       pattern = adaptor[adaptor.length - i ... adaptor.length]
 123       mis     = (pattern.length * misp * 0.01).round
 124       ins     = (pattern.length * insp * 0.01).round
 125       del     = (pattern.length * delp * 0.01).round
 126
 127       if i == adaptor.length
 128         patterns << "#{pattern}[#{mis},#{ins},#{del}]"
 129       else
 130         patterns << "^ #{pattern}[#{mis},#{ins},#{del}]"
 131       end
 132     end
 133
 134     compile(patterns)
 135   end
 136
 137   def self.create_right(adaptor, misp, insp, delp, len)
 138     patterns = []
 139
 140     (adaptor.length).downto(len) do |i|
 141       pattern = adaptor[0 ... i]
 142       mis     = (pattern.length * misp * 0.01).round
 143       ins     = (pattern.length * insp * 0.01).round
 144       del     = (pattern.length * delp * 0.01).round
 145
 146       if i == adaptor.length
 147         patterns << "#{pattern}[#{mis},#{ins},#{del}]"
 148       else
 149         patterns << "#{pattern}[#{mis},#{ins},#{del}] $"
 150       end
 151     end
 152
 153     compile(patterns)
 154   end
 155
 156   private
 157
 158   def self.compile(patterns)
 159     new_patterns = []
 160
 161     while patterns.size > 1
 162       new_patterns = patterns[0 ... -2]
 163       new_patterns << "( #{patterns[-2 .. -1].join(' | ')} )"
 164
 165       patterns = new_patterns
 166     end
 167
 168     patterns.first
 169   end
 170 end
 171
 172 casts = []
 173 casts << {:long=>'forward',       :short=>'f', :type=>'string', :mandatory=>false, :default=>nil,        :allowed=>nil, :disallowed=>nil}
 174 casts << {:long=>'forward_rc',    :short=>'F', :type=>'string', :mandatory=>false, :default=>nil,        :allowed=>nil, :disallowed=>nil}
 175 casts << {:long=>'reverse',       :short=>'r', :type=>'string', :mandatory=>false, :default=>nil,        :allowed=>nil, :disallowed=>nil}
 176 casts << {:long=>'reverse_rc',    :short=>'R', :type=>'string', :mandatory=>false, :default=>nil,        :allowed=>nil, :disallowed=>nil}
 177 casts << {:long=>'len_forward',   :short=>'l', :type=>'uint',   :mandatory=>false, :default=>nil,        :allowed=>nil, :disallowed=>'0'}
 178 casts << {:long=>'len_reverse',   :short=>'L', :type=>'uint',   :mandatory=>false, :default=>nil,        :allowed=>nil, :disallowed=>'0'}
 179 casts << {:long=>'mismatches',    :short=>'m', :type=>'uint',   :mandatory=>false, :default=>10,         :allowed=>nil, :disallowed=>nil}
 180 casts << {:long=>'insertions',    :short=>'i', :type=>'uint',   :mandatory=>false, :default=>5,          :allowed=>nil, :disallowed=>nil}
 181 casts << {:long=>'deletions',     :short=>'d', :type=>'uint',   :mandatory=>false, :default=>5,          :allowed=>nil, :disallowed=>nil}
 182 casts << {:long=>'cpus',          :short=>'c', :type=>'uint',   :mandatory=>false, :default=>1,          :allowed=>nil, :disallowed=>'0'}
 183 casts << {:long=>'bases_per_file',:short=>'b', :type=>'uint',   :mandatory=>false, :default=>10_000_000, :allowed=>nil, :disallowed=>'0'}
 184
 185 options = Biopieces.options_parse(ARGV, casts)
 186
 187 if options[:forward_rc]
 188   options[:forward] = Seq.new("test", options[:forward_rc], 'dna').revcomp.seq
 189 end
 190
 191 if options[:reverse_rc]
 192   options[:reverse] = Seq.new("test", options[:reverse_rc], 'dna').revcomp.seq
 193 end
 194
 195 raise ArgumentError, "no adaptor specified" unless options[:forward] or options[:reverse]
 196
 197 if options[:forward]
 198   options[:len_forward] = options[:forward].length unless options[:len_forward]
 199
 200   if options[:len_forward] > options[:forward].length
 201     raise ArgumentError, "len_forward > forward adaptor (#{options[:len_forward]} > #{options[:forward].length})"
 202   end
 203 end
 204
 205 if options[:reverse]
 206   options[:len_reverse] = options[:reverse].length unless options[:len_reverse]
 207
 208   if options[:len_reverse] > options[:reverse].length
 209     raise ArgumentError, "len_reverse > reverse adaptor (#{options[:len_reverse]} > #{options[:reverse].length})"
 210   end
 211 end
 212
 213 tmpdir       = Biopieces.mktmpdir
 214 file_records = File.join(tmpdir, "data.stream")
 215
 216 count_file = 0
 217 count_seq  = 0
 218 bases      = 0
 219
 220 Biopieces.open(options[:stream_in], file_records) do |input, output|
 221   file_fasta = File.join(tmpdir, "#{count_file}.fna")
 222   out_fa     = Fasta.open(file_fasta, 'w')
 223
 224   input.each do |record|
 225     output.puts record
 226
 227     if record[:SEQ] and record[:SEQ].length > 0
 228       record[:SEQ_NAME] = count_seq
 229
 230       seq = Seq.new_bp(record)
 231
 232       out_fa.puts seq.to_fasta
 233
 234       count_seq += 1;
 235       bases      += record[:SEQ].length
 236
 237       if bases > options[:bases_per_file]
 238         out_fa.close
 239         bases       = 0
 240         count_file += 1
 241         file_fasta  = File.join(tmpdir, "#{count_file}.fna")
 242         out_fa      = Fasta.open(file_fasta, 'w')
 243       end
 244     end
 245   end
 246
 247   out_fa.close if out_fa.respond_to? :close
 248 end
 249
 250 patscan       = PatScan.new(tmpdir, options[:cpus])
 251 matches_left  = {}
 252 matches_right = {}
 253
 254 if options[:forward]
 255   pat_left = Pattern.create_left(
 256     options[:forward],
 257     options[:mismatches],
 258     options[:insertions],
 259     options[:deletions],
 260     options[:len_forward])
 261
 262   patscan.run(pat_left)
 263   matches_left = patscan.parse_results
 264 end
 265
 266 if options[:reverse]
 267   patscan.clean
 268   pat_right = Pattern.create_right(
 269     options[:reverse],
 270     options[:mismatches],
 271     options[:insertions],
 272     options[:deletions],
 273     options[:len_reverse])
 274
 275   patscan.run(pat_right)
 276   matches_right = patscan.parse_results
 277 end
 278
 279 count_seq = 0
 280
 281 Biopieces.open(file_records, options[:stream_out]) do |input, output|
 282   input.each_record do |record|
 283     if record[:SEQ]
 284       if matches_left[count_seq]
 285         record[:ADAPTOR_POS_LEFT] = matches_left[count_seq].first
 286         record[:ADAPTOR_LEN_LEFT] = matches_left[count_seq].last
 287         record[:ADAPTOR_PAT_LEFT] = record[:SEQ][record[:ADAPTOR_POS_LEFT] ... record[:ADAPTOR_POS_LEFT] + record[:ADAPTOR_LEN_LEFT]]
 288       end
 289
 290       if matches_right[count_seq]
 291         record[:ADAPTOR_POS_RIGHT] = matches_right[count_seq].first
 292         record[:ADAPTOR_LEN_RIGHT] = matches_right[count_seq].last
 293         record[:ADAPTOR_PAT_RIGHT] = record[:SEQ][record[:ADAPTOR_POS_RIGHT] ... record[:ADAPTOR_POS_RIGHT] + record[:ADAPTOR_LEN_RIGHT]]
 294       end
 295
 296       count_seq += 1;
 297     end
 298
 299     output.puts record
 300   end
 301 end
 302
 303
 304 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 305
 306
 307 __END__