bp_bin/uclust_seq

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2010 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 #
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31
  32 require 'biopieces'
  33 require 'fasta'
  34
  35 class Uclust
  36   include Enumerable
  37
  38   def initialize(infile, outfile)
  39     @infile  = infile
  40     @outfile = outfile
  41   end
  42
  43   # Method that calls Uclusts sorting for sorting a FASTA file
  44   # according to decending sequence length.
  45   def sort(options)
  46     command = "nice -n 19"
  47     command << " uclust --sort #{@infile} --output #{@infile}.sort"
  48     command << " > /dev/null 2>&1" unless options[:verbose]
  49     system(command)
  50     raise "Failed sorting file: #{@infile}" unless $?.success?
  51
  52     File.rename "#{@infile}.sort", @infile
  53   end
  54
  55   def ublast(options)
  56     # uclust --ublast query_seqs.fasta --db database.fasta --blast6out filename --evalue E
  57     options[:e_val] = 10 if options[:e_val].is_nil?
  58     command = "nice -n 19"
  59     command << " uclust --ublast #{@infile} --db #{options[:database]} --uc #{@outfile} --evalue #{options[:e_val]}"
  60     command << " > /dev/null 2>&1" unless options[:verbose]
  61     system(command)
  62     raise "Command failed: #{command}" unless $?.success?
  63   end
  64
  65   def usearch(options)
  66     # uclust --query query.fasta --db db.fasta --uc results.uc --id 0.90 [--evalue E]
  67     command = "nice -n 19"
  68     command << " uclust --query #{@infile} --db #{options[:database]} --uc #{@outfile} --id #{options[:identity]}"
  69     command << " --evalue #{options[:e_val]}" if options.has_key? :e_val
  70     command << " > /dev/null 2>&1" unless options[:verbose]
  71     system(command)
  72     raise "Command failed: #{command}" unless $?.success?
  73   end
  74
  75   def uclust(options)
  76     # uclust --input seqs_sorted.fasta --uc results.uc --id 0.90
  77     command = "nice -n 19"
  78     command << " uclust --input #{@infile} --uc #{@outfile} --id #{options[:identity]}"
  79     command << " > /dev/null 2>&1" unless options[:verbose]
  80     system(command)
  81     raise "Command failed: #{command}" unless $?.success?
  82   end
  83
  84   def usearch_uclust(options)
  85     # uclust --input seqs_sorted.fasta --lib db.fasta --uc results.uc --id 0.90
  86     command = "nice -n 19"
  87     command << " uclust --input #{@infile} --lib #{options[:database]} --uc #{@outfile} --id #{options[:identity]}"
  88     command << " --lib #{options[:database]}" if options.has_key? :database
  89     command << " > /dev/null 2>&1" unless options[:verbose]
  90     system(command)
  91     raise "Command failed: #{command}" unless $?.success?
  92   end
  93
  94   def each
  95     record = {}
  96
  97     File.open(@outfile, mode="r") do |ios|
  98       ios.each_line do |line|
  99         if line !~ /^#/
 100           fields = line.split("\t")
 101
 102           record[:REC_TYPE] = "UCLUST"
 103           record[:TYPE]     = fields[0]
 104           record[:CLUSTER]  = fields[1]
 105           record[:SEQ_LEN]  = fields[2]
 106           record[:IDENT]    = fields[3]
 107           record[:STRAND]   = fields[4]
 108           record[:Q_BEG]    = fields[5]
 109           record[:S_BEG]    = fields[6]
 110           record[:CIGAR]    = fields[7]
 111           record[:Q_ID]     = fields[8]
 112           record[:S_ID]     = fields[9]
 113
 114           yield record
 115         end
 116       end
 117     end
 118
 119     self # conventionally
 120   end
 121 end
 122
 123 ok_methods = "ublast,usearch,uclust,usearch_uclust"
 124
 125 casts = []
 126 casts << {:long=>'no_sort',  :short=>'n', :type=>'flag',   :mandatory=>false, :default=>nil,      :allowed=>nil,        :disallowed=>nil}
 127 casts << {:long=>'method',   :short=>'m', :type=>'string', :mandatory=>true,  :default=>"uclust", :allowed=>ok_methods, :disallowed=>nil}
 128 casts << {:long=>'database', :short=>'d', :type=>'file!',  :mandatory=>false, :default=>nil,      :allowed=>nil,        :disallowed=>nil}
 129 casts << {:long=>'identity', :short=>'i', :type=>'float',  :mandatory=>true,  :default=>0.9,      :allowed=>nil,        :disallowed=>nil}
 130 casts << {:long=>'e_val',    :short=>'e', :type=>'float',  :mandatory=>false, :default=>nil,      :allowed=>nil,        :disallowed=>nil}
 131
 132 bp = Biopieces.new
 133
 134 options = bp.parse(ARGV, casts)
 135
 136 tmpdir  = bp.mktmpdir
 137 infile  = "#{tmpdir}/in.fna"
 138 outfile = "#{tmpdir}/out.uc"
 139
 140 Fasta.open(infile, mode="w") do |fasta_io|
 141   bp.each_record do |record|
 142     bp.puts record
 143     fasta_io.puts record
 144   end
 145 end
 146
 147 uclust = Uclust.new(infile, outfile)
 148 uclust.sort(options) unless options[:no_sort]
 149
 150 case options[:method].to_s
 151 when "ublast"         then uclust.ublast(options)
 152 when "usearch"        then uclust.usearch(options)
 153 when "uclust"         then uclust.uclust(options)
 154 when "usearch_uclust" then uclust.usearch_uclust(options)
 155 else raise "Unknown method: #{options[:method]}"
 156 end
 157
 158 uclust.each do |record|
 159   bp.puts record
 160 end
 161
 162
 163 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 164
 165
 166 __END__