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