# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-require 'biopieces'
-require 'seq'
-
-VERBOSE = true
-
-def disambiguate(adaptor)
- adaptor_disamb = adaptor.dup
- adaptor_disamb.gsub!('U', 'T')
- adaptor_disamb.gsub!('R', '[AG]')
- adaptor_disamb.gsub!('Y', '[CT]')
- adaptor_disamb.gsub!('S', '[GC]')
- adaptor_disamb.gsub!('W', '[AT]')
- adaptor_disamb.gsub!('M', '[AC]')
- adaptor_disamb.gsub!('K', '[GT]')
- adaptor_disamb.gsub!('V', '[ACG]')
- adaptor_disamb.gsub!('H', '[ACT]')
- adaptor_disamb.gsub!('D', '[AGT]')
- adaptor_disamb.gsub!('B', '[CGT]')
- adaptor_disamb.gsub!('N', '.')
- adaptor_disamb
-end
+require 'pp'
+require 'maasha/biopieces'
+require 'maasha/fasta'
+
+# Error class for PatScan errors.
+class PatScanError < StandardError; end;
+
+class PatScan
+ def initialize(options, tmpdir, file_pattern, cpus)
+ @options = options
+ @tmpdir = tmpdir
+ @file_pattern = file_pattern
+ @cpus = cpus
+ @files_fasta = Dir.glob(File.join(@tmpdir, "*.fna"))
+
+ pat = Pattern.new(@options)
+ pat.write(@file_pattern)
+ end
-class Seq
- # Method that finds an adaptor or part thereof in the sequence of a Seq object.
- # Returns a Match object if the adaptor was found otherwise nil. The ed_percent
- # indicates the maximum edit distance allowed in all possible overlaps.
- def adaptor_find(adaptor, adaptor_disamb, pos = 0, ed_percent = 0)
- raise SeqError, "Edit distance percent out of range #{ed_percent}" unless (0 .. 100).include? ed_percent
+ def run
+ child_count = 0
+ ch_mutex = Mutex.new
+ threads = []
- if pos < 0
- pos = self.length + pos # pos offset from the right end
- end
+ @files_fasta.each do |file|
+ Thread.pass while child_count >= @cpus
+ ch_mutex.synchronize { child_count += 1 }
- if pos < self.length
- if match = adaptor_find_simple(adaptor_disamb, pos)
- return match
- elsif match = adaptor_find_complex(adaptor, pos, ed_percent)
- return match
- elsif match = adaptor_partial_find_complex(adaptor, pos, ed_percent)
- return match
+ threads << Thread.new do
+ command = command_compile(file)
+ system(command)
+ raise PatScanError, "Command failed: #{command}" unless $?.success?
+ ch_mutex.synchronize { child_count -= 1 }
end
end
+
+ threads.each { |t| t.join }
end
- private
+ def command_compile(file)
+ commands = []
+ commands << "nice -n 19"
+ commands << "scan_for_matches"
+ commands << @file_pattern
+ commands << "< #{file}"
+ commands << "> #{file}.out"
+ command = commands.join(" ")
+ end
- # Method to find an adaptor in a sequence taking into account ambiguity
- # codes, but not considering mismatches, insertions, and deletions.
- def adaptor_find_simple(adaptor, pos)
- self.seq.upcase.match(adaptor, pos) do |m|
- return Match.new($`.length, m, m.to_s.length, 0, 0, 0, m.to_s.length)
+ def parse_results
+ files_result = Dir.glob(File.join(@tmpdir, "*.out"))
+
+ matches = {}
+
+ files_result.each do |file|
+ Fasta.open(file, 'r') do |ios|
+ ios.each do |entry|
+ if entry.seq_name =~ /^(\d+):\[(\d+),(\d+)\]$/
+ name = $1.to_i
+ start = $2.to_i - 1
+ stop = $3.to_i - 1
+ matches[name] = [start, stop - start + 1] unless matches.has_key? name
+ else
+ raise "Failed to parse sequence name: #{entry.seq_name}"
+ end
+ end
+ end
end
- end
- # Method to find an adaptor in a sequence taking into account ambiguity
- # codes, mismatches, insertions, and deletions.
- def adaptor_find_complex(adaptor, pos, ed_percent)
- ed_max = (adaptor.length * ed_percent * 0.01).round
+ matches
+ end
+end
- match = self.match(adaptor, pos, ed_max)
+# Error class for Pattern errors.
+class PatternError < StandardError; end;
- match
+class Pattern
+ def initialize(options)
+ @options = options
+ @patterns = []
+ @patterns << pattern_internal
+ @patterns += patterns_end if @options[:partial]
end
- # Method to find part of an adaptor at the right end of a sequence taking
- # into account ambiguity codes, mismatches, insertions, and deletions.
- def adaptor_partial_find_complex(adaptor, pos, ed_percent)
- if pos > self.length - adaptor.length
- adaptor = adaptor[0 ... self.length - pos]
- else
- adaptor = adaptor[0 ... adaptor.length - 1]
+ def to_i
+ new_patterns = []
- pos = self.length - adaptor.length
+ while @patterns.size > 1
+ new_patterns = @patterns[0 ... -2]
+ new_patterns << "( #{@patterns[-2 .. -1].join(' | ')} )"
+
+ @patterns = new_patterns
end
- #puts self.seq if VERBOSE
+ @patterns.first
+ end
+
+ def write(file)
+ File.open(file, 'w') do |ios|
+ ios.puts self.to_i
+ end
+ end
- while adaptor.length > 0
- #puts (" " * pos) + adaptor if VERBOSE
+ private
- ed_max = (adaptor.length * ed_percent * 0.01).round
+ def pattern_internal
+ pattern = @options[:adaptor]
+ mis = mis_count(pattern)
+ ins = ins_count(pattern)
+ del = del_count(pattern)
- if ed_max == 0
- self.seq.upcase.match(adaptor, pos) do |m|
- return Match.new($`.length, m, m.to_s.length, 0, 0, 0, m.to_s.length)
- end
- else
- self.scan(adaptor, pos, ed_max).each do |match|
- return match
- end
- end
+ "#{pattern}[#{mis},#{ins},#{del}]"
+ end
- adaptor = adaptor[0 ... -1]
+ def patterns_end
+ patterns = []
+ adaptor = @options[:adaptor]
- pos += 1
+ raise PatternError, "len > adaptor length: #{@options[:len]} > #{adaptor.length - 1}" if @options[:len] > adaptor.length - 1
+
+ (adaptor.length - 1).downto(@options[:len]) do |i|
+ pattern = adaptor[0 ... i]
+ mis = mis_count(pattern)
+ ins = ins_count(pattern)
+ del = del_count(pattern)
+ patterns << "#{pattern}[#{mis},#{ins},#{del}] $"
end
+
+ patterns
+ end
+
+ def mis_count(pattern)
+ (pattern.length * @options[:mismatches] * 0.01).round
+ end
+
+ def ins_count(pattern)
+ (pattern.length * @options[:insertions] * 0.01).round
+ end
+
+ def del_count(pattern)
+ (pattern.length * @options[:deletions] * 0.01).round
end
end
casts = []
-casts << {:long=>'adaptor', :short=>'r', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
-casts << {:long=>'edit_distance', :short=>'e', :type=>'uint', :mandatory=>false, :default=>20, :allowed=>nil, :disallowed=>nil}
-casts << {:long=>'pos', :short=>'p', :type=>'int', :mandatory=>false, :default=>1, :allowed=>nil, :disallowed=>"0"}
-casts << {:long=>'cache', :short=>'c', :type=>'flag', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'adaptor', :short=>'a', :type=>'string', :mandatory=>true, :default=>nil, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'partial', :short=>'p', :type=>'flag', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'len', :short=>'l', :type=>'uint', :mandatory=>false, :default=>10, :allowed=>nil, :disallowed=>'0'}
+casts << {:long=>'mismatches', :short=>'m', :type=>'uint', :mandatory=>false, :default=>10, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'insertions', :short=>'i', :type=>'uint', :mandatory=>false, :default=>5, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'deletions', :short=>'d', :type=>'uint', :mandatory=>false, :default=>5, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'cpus', :short=>'c', :type=>'uint', :mandatory=>false, :default=>1, :allowed=>nil, :disallowed=>'0'}
-bp = Biopieces.new
+BASE_PER_FILE = 10_000_000
-options = bp.parse(ARGV, casts)
+options = Biopieces.options_parse(ARGV, casts)
-adaptor = options[:adaptor].to_s.upcase
-adaptor_disamb = disambiguate(adaptor)
+tmpdir = Biopieces.mktmpdir
+file_records = File.join(tmpdir, "data.stream")
+file_pattern = File.join(tmpdir, "pattern.txt")
-pos = options[:pos]
-pos -= 1 if pos > 0 # pos was 1-based
+number_file = 0
+number_seq = 0
+bases = 0
-cache = {}
+Biopieces.open(options[:stream_in], file_records) do |input, output|
+ file_fasta = File.join(tmpdir, "#{number_file}.fna")
+ out_fa = Fasta.open(file_fasta, 'w')
-bp.each_record do |record|
- if record.has_key? :SEQ
- entry = Seq.new(record[:SEQ_NAME], record[:SEQ], "dna", record[:SCORES])
+ input.each do |record|
+ output.puts record
- if cache[entry.seq.upcase.to_sym] and options[:cache]
- match = cache[entry.seq.upcase]
- else
- match = entry.adaptor_find(adaptor, adaptor_disamb, pos, options[:edit_distance])
+ if record.has_key? :SEQ
+ record[:SEQ_NAME] = number_seq
- cache[entry.seq.upcase.to_sym] = match if match and options[:cache]
- end
+ seq = Seq.new_bp(record)
+
+ out_fa.puts seq.to_fasta
+
+ number_seq += 1;
+ bases += record[:SEQ].length
- if match
- record[:ADAPTOR_POS] = match.pos
- record[:ADAPTOR_LEN] = match.length
- record[:ADAPTOR_MATCH] = match.match
+ if bases > BASE_PER_FILE
+ out_fa.close
+ bases = 0
+ number_file += 1
+ file_fasta = File.join(tmpdir, "#{number_file}.fna")
+ out_fa = Fasta.open(file_fasta, 'w')
+ end
end
end
- bp.puts record
+ out_fa.close if out_fa.respond_to? :close
+end
+
+patscan = PatScan.new(options, tmpdir, file_pattern, options[:cpus])
+patscan.run
+matches = patscan.parse_results
+
+number_seq = 0
+
+Biopieces.open(file_records, options[:stream_out]) do |input, output|
+ input.each_record do |record|
+ if record.has_key? :SEQ
+ if matches.has_key? number_seq
+ record[:ADAPTOR_POS] = matches[number_seq].first
+ record[:ADAPTOR_LEN] = matches[number_seq].last
+ end
+
+ number_seq += 1;
+ end
+
+ output.puts record
+ end
end
projects / biopieces.git / blobdiff