added requirement of RubyInline

[biopieces.git] / bp_bin / uclust_seq

diff --git a/bp_bin/uclust_seq b/bp_bin/uclust_seq
index 4523ab05b09bf69a6919b9ac9009bed60bbd0bc1..b20a29c672673ce690870e7f1e31d54df24e30b0 100755 (executable)
--- a/bp_bin/uclust_seq
+++ b/bp_bin/uclust_seq
@@ -1,6 +1,6 @@
 #!/usr/bin/env ruby
 
-# Copyright (C) 2007-2010 Martin A. Hansen.
+# Copyright (C) 2007-2011 Martin A. Hansen.
 
 # This program is free software; you can redistribute it and/or
 # modify it under the terms of the GNU General Public License
@@ -28,10 +28,11 @@
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-
 require 'maasha/biopieces'
 require 'maasha/fasta'
 
+SORT_LIMIT = 2_000_000_000   # use mergesort for files biggern than 2Gb.
+
 class Uclust
   include Enumerable
 
@@ -42,46 +43,41 @@ class Uclust
     @command = []
   end
 
-  # Method that calls Uclusts sorting for sorting a FASTA file
+  # Method that calls Usearch sorting for sorting a FASTA file
   # according to decending sequence length.
   def sort
-    @command << "uclust --sort #{@infile} --output #{@infile}.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
 
-  def ublast
-    # uclust --ublast query_seqs.fasta --db database.fasta --blast6out filename --evalue E
-    @options[:e_val] = 10 unless @options[:e_val]
-    @command << "uclust --ublast #{@infile} --db #{@options[:database]} --blast6out #{@outfile} --evalue #{@options[:e_val]}"
+       # Method to execute clustering de novo.
+  def cluster
+    @command << "usearch --cluster #{@infile} --uc #{@outfile} --id #{@options[:identity]}"
 
                execute
   end
 
   # Method to execute database search.
   def usearch
-    # uclust --query query.fasta --db db.fasta --uc results.uc --id 0.90 [--evalue E]
-    @command << "uclust --query #{@infile} --db #{@options[:database]} --uc #{@outfile} --id #{@options[:identity]}"
+    # usearch --query query.fasta --db db.fasta --uc results.uc --id 0.90 [--evalue E]
+    @command << "usearch --query #{@infile} --db #{@options[:database]} --uc #{@outfile} --id #{@options[:identity]}"
     @command << "--evalue #{@options[:e_val]}" if @options.has_key? :e_val
 
                execute
   end
 
-       # Method to execute clustering de novo.
-  def uclust
-    # uclust --input seqs_sorted.fasta --uc results.uc --id 0.90
-    @command << "uclust --input #{@infile} --uc #{@outfile} --id #{@options[:identity]}"
-
-               execute
-  end
-
        # Method to execute clustering to database plus de novo if not matched.
   def usearch_uclust
-    # uclust --input seqs_sorted.fasta --lib db.fasta --uc results.uc --id 0.90
-    @command << "uclust --input #{@infile} --lib #{@options[:database]} --uc #{@outfile} --id #{@options[:identity]}"
-    @command << "--lib #{@options[:database]}" if @options.has_key? :database
+    # usearch --cluster seqs_sorted.fasta --db db.fasta --uc results.uc --id 0.90
+    @command << "usearch --cluster #{@infile} --db #{@options[:database]} --uc #{@outfile} --id #{@options[:identity]}"
 
                execute
   end
@@ -98,12 +94,13 @@ class Uclust
 
           record[:REC_TYPE] = "UCLUST"
           record[:TYPE]     = fields[0]
-          record[:CLUSTER]  = fields[1]
-          record[:SEQ_LEN]  = fields[2]
-          record[:IDENT]    = fields[3]
+          record[:CLUSTER]  = fields[1].to_i
+          record[:SEQ_LEN]  = fields[2].to_i
+          record[:IDENT]    = fields[3].to_f
           record[:STRAND]   = fields[4]
-          record[:Q_BEG]    = fields[5]
-          record[:S_BEG]    = fields[6]
+          record[:Q_BEG]    = fields[5].to_i
+          record[:S_BEG]    = fields[6].to_i
+          record[:S_END]    = fields[6].to_i + fields[2].to_i
           record[:CIGAR]    = fields[7]
           record[:Q_ID]     = fields[8]
           record[:S_ID]     = fields[9]
@@ -132,7 +129,7 @@ class Uclust
        end
 end
 
-ok_methods = "ublast,usearch,uclust,usearch_uclust"
+ok_methods = "uclust,usearch,usearch_uclust"
 
 casts = []
 casts << {:long=>'no_sort',  :short=>'n', :type=>'flag',   :mandatory=>false, :default=>nil,      :allowed=>nil,        :disallowed=>nil}
@@ -144,6 +141,10 @@ casts << {:long=>'e_val',    :short=>'e', :type=>'float',  :mandatory=>false, :d
 
 options = Biopieces.options_parse(ARGV, casts)
 
+# At high identities, around 96% and above, compressed indexes are often more sensitive, faster 
+# and use less RAM. Compressed indexes are disabled by default, so I generally recommend that 
+# you specify the --slots and --w options when clustering at high identities.
+
 tmpdir  = Biopieces.mktmpdir
 infile  = File.join(tmpdir, "in.fna")
 outfile = File.join(tmpdir, "out.uc")
@@ -152,19 +153,20 @@ Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
   Fasta.open(infile, mode="w") do |fasta_io|
     input.each_record do |record|
       output.puts record
-      fasta_io.puts record
+
+      if record.has_key? :SEQ_NAME and record.has_key? :SEQ
+        fasta_io.puts Seq.new_bp(record).to_fasta
+      end
     end
   end
 
   uclust = Uclust.new(infile, outfile, options)
-  uclust.sort unless options[:no_sort]
+  uclust.sort unless options[:no_sort] or options[:method] == "usearch"
 
   case options[:method].to_s
-  when "ublast"         then uclust.ublast
+  when "uclust"         then uclust.cluster
   when "usearch"        then uclust.usearch
-  when "uclust"         then uclust.uclust
   when "usearch_uclust" then uclust.usearch_uclust
-  else raise "Unknown method: #{options[:method]}"
   end
 
   uclust.each do |record|