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