]> git.donarmstrong.com Git - biopieces.git/blobdiff - bp_bin/uclust_seq
added missing cpus switch to uclust_seq
[biopieces.git] / bp_bin / uclust_seq
index 0eaee60bc66ba6fd54b15752d7c394d89fa28483..5a6d0ef901a1f9cb298b324106c568c721a50b1d 100755 (executable)
 
 require 'maasha/biopieces'
 require 'maasha/fasta'
-
-SORT_LIMIT = 2_000_000_000   # use mergesort for files biggern than 2Gb.
-
-class Uclust
-  include Enumerable
-
-  def initialize(infile, outfile, options)
-    @infile  = infile
-    @outfile = outfile
-    @options = options
-    @command = []
-  end
-
-  # Method that calls Usearch sorting for sorting a FASTA file
-  # according to decending sequence length.
-  def sort
-    # usearch -sort seqs.fasta -output seqs.sorted.fasta
-    if File.size(@infile) < SORT_LIMIT
-      @command << "usearch --sort #{@infile} --output #{@infile}.sort"
-    else
-      @command << "usearch --mergesort #{@infile} --output #{@infile}.sort"
-    end
-
-               execute
-
-    File.rename "#{@infile}.sort", @infile
-  end
-
-       # Method to execute clustering de novo.
-  def cluster
-    @command << "usearch --cluster #{@infile} --uc #{@outfile} --id #{@options[:identity]}"
-
-               execute
-  end
-
-       # Method to parse a Uclust .uc file and for each line of data
-       # yield a Biopiece record.
-  def each
-    record = {}
-
-    File.open(@outfile, mode="r") do |ios|
-      ios.each_line do |line|
-        if line !~ /^#/
-          fields = line.chomp.split("\t")
-
-          next if fields[0] == 'C'
-
-          record[:TYPE]     = fields[0]
-          record[:CLUSTER]  = fields[1].to_i
-          record[:IDENT]    = fields[3].to_f
-          record[:Q_ID]     = fields[8]
-
-          yield record
-        end
-      end
-    end
-
-    self # conventionally
-  end
-
-  private
-
-       # Method to execute a command using a system() call.
-       # The command is composed of bits from the @command variable.
-       def execute
-               @command.unshift "nice -n 19"
-    @command << "--rev" if @options[:comp]
-               @command << "> /dev/null 2>&1" unless @options[:verbose]
-               command = @command.join(" ")
-    system(command)
-    raise "Command failed: #{command}" unless $?.success?
-
-               @command = []
-       end
-end
+require 'maasha/usearch'
 
 casts = []
-casts << {:long=>'no_sort',  :short=>'n', :type=>'flag',   :mandatory=>false, :default=>nil,      :allowed=>nil,        :disallowed=>nil}
-casts << {:long=>'comp',     :short=>'c', :type=>'flag',   :mandatory=>false, :default=>nil,      :allowed=>nil,        :disallowed=>nil}
-casts << {:long=>'identity', :short=>'i', :type=>'float',  :mandatory=>true,  :default=>0.9,      :allowed=>nil,        :disallowed=>nil}
+casts << {long: 'no_sort',  short: 'n', type: 'flag',  mandatory: false, default: nil, allowed: nil, disallowed: nil}
+casts << {long: 'comp',     short: 'c', type: 'flag',  mandatory: false, default: nil, allowed: nil, disallowed: nil}
+casts << {long: 'identity', short: 'i', type: 'float', mandatory: true,  default: 0.9, allowed: nil, disallowed: nil}
+casts << {long: 'cpus',     short: 'C', type: 'uint',  mandatory: false, default: 1,   allowed: nil, disallowed: "0"}
 
 options = Biopieces.options_parse(ARGV, casts)
 
@@ -119,35 +46,36 @@ file_fasta   = File.join(tmpdir, "in.fna")
 file_uclust  = File.join(tmpdir, "out.uc")
 
 Biopieces.open(options[:stream_in], file_records) do |input, output|
-  Fasta.open(file_fasta, mode="w") do |fasta_io|
+  Fasta.open(file_fasta, 'w') do |fasta_io|
     input.each_record do |record|
       output.puts record
 
-      if record.has_key? :SEQ_NAME and record.has_key? :SEQ
+      if record[:SEQ_NAME] and record[:SEQ]
         fasta_io.puts Seq.new_bp(record).to_fasta
       end
     end
   end
 end
 
-uc = Uclust.new(file_fasta, file_uclust, options)
-uc.sort unless options[:no_sort]
-uc.cluster
+us = Usearch.new(file_fasta, file_uclust, options)
+
+us.sortbylength unless options[:no_sort]
+us.cluster_smallmem
 
 hash = {}
 
-uc.each do |record|
-  hash[record[:Q_ID].to_sym] = record.dup
+us.each_cluster do |cluster|
+  hash[cluster[:Q_ID].to_sym] = cluster.dup
 end
 
 Biopieces.open(file_records, options[:stream_out]) do |input, output|
   input.each_record do |record|
-    if record.has_key? :SEQ_NAME and record.has_key? :SEQ
-      if hash.has_key? record[:SEQ_NAME].to_sym
-        uc = hash[record[:SEQ_NAME].to_sym]
-        record[:CLUSTER] = uc[:CLUSTER].to_i
-        record[:IDENT]   = uc[:IDENT].to_i
-        record[:IDENT] = '*' if uc[:TYPE] == 'S'
+    if record[:SEQ_NAME] and record[:SEQ]
+      if hash[record[:SEQ_NAME].to_sym]
+        us = hash[record[:SEQ_NAME].to_sym]
+        record[:CLUSTER] = us[:CLUSTER].to_i
+        record[:IDENT]   = us[:IDENT].to_i
+        record[:IDENT] = '*' if us[:TYPE] == 'S'
       end
     end