X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bp_bin%2Fusearch_seq;h=946659f0c165652d32f2c39259ef8074e9484be8;hb=bdcd4e373c9a3361781e5f171d2d60166c46dc3d;hp=e563893f6eaa2ad4a02b959282c04ed7bc89bc8b;hpb=a1dea64932bbda4eb2a90de9d3e5d353ebeed295;p=biopieces.git

diff --git a/bp_bin/usearch_seq b/bp_bin/usearch_seq
index e563893..946659f 100755
--- a/bp_bin/usearch_seq
+++ b/bp_bin/usearch_seq
@@ -1,6 +1,6 @@
-#!/usr/bin/env perl
+#!/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
@@ -18,205 +18,64 @@
 
 # http://www.gnu.org/copyleft/gpl.html
 
-
-# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-# USEARCH sequences in the stream against a specified database or genome.
-
-# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-
-use warnings;
-use strict;
-use Data::Dumper;
-use Maasha::Biopieces;
-use Maasha::Common;
-use Maasha::Seq;
-use Maasha::Fasta;
-
-
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
+# This program is part of the Biopieces framework (www.biopieces.org).
 
-my ( $TYPE_HASH, $options, $in, $out, $tmp_dir, $tmp_in, $tmp_out, $q_type, $record, $entry, $fh_in, $fh_out );
-
-$TYPE_HASH = {
-    'L' => 'LibSeed',
-    'S' => 'NewSeed',
-    'H' => 'Hit',
-    'R' => 'Reject',
-    'D' => 'LibCluster',
-    'C' => 'NewCluster',
-    'N' => 'NoHit'
-};
-
-$options = Maasha::Biopieces::parse_options(
-    [
-        { long => 'database', short => 'd', type => 'file',   mandatory => 'no', default => undef,  allowed => undef, disallowed => undef },
-        { long => 'genome',   short => 'g', type => 'genome', mandatory => 'no', default => undef,  allowed => undef, disallowed => undef },
-        { long => 'id',       short => 'i', type => 'float',  mandatory => 'no', default => '0.90', allowed => undef, disallowed => undef },
-    ]   
-);
-
-Maasha::Common::error( qq(both --database and --genome specified) ) if     $options->{ "genome" } and     $options->{ "database" };
-Maasha::Common::error( qq(no --database or --genome specified) )    if not $options->{ "genome" } and not $options->{ "database" };
-
-$in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
-$out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
-
-$options->{ "database" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna" if $options->{ 'genome' };
-
-$tmp_dir = Maasha::Biopieces::get_tmpdir();
-$tmp_in  = "$tmp_dir/usearch_query.seq";
-$tmp_out = "$tmp_dir/usearch.uc";
-
-$fh_out = Maasha::Filesys::file_write_open( $tmp_in );
-
-while ( $record = Maasha::Biopieces::get_record( $in ) ) 
-{
-    if ( $entry = Maasha::Fasta::biopiece2fasta( $record ) )
-    {
-        $q_type = Maasha::Seq::seq_guess_type( $record->{ 'SEQ' } ) if not $q_type;
-
-        Maasha::Fasta::put_entry( $entry, $fh_out );
-    }
-
-    Maasha::Biopieces::put_record( $record, $out );
-}
-
-close $fh_out;
-
-
-if ( $options->{ 'verbose' } )
-{
-    Maasha::Common::run(
-        "uclust",
-        join( " ",
-            "--input $tmp_in",
-            "--lib $options->{ 'database' }",
-            "--libonly",
-            "--uc $tmp_out",
-            "--id $options->{ 'id' }",
-            "--rev",
-        ),
-        1
-    );
-}
-else
-{
-    Maasha::Common::run(
-        "uclust",
-        join( " ",
-            "--input $tmp_in",
-            "--lib $options->{ 'database' }",
-            "--libonly",
-            "--uc $tmp_out",
-            "--id $options->{ 'id' }",
-            "--rev",
-            "> /dev/null 2>&1"
-        ),
-        1
-    );
-}
-
-unlink $tmp_in;
-
-$fh_in = Maasha::Filesys::file_read_open( $tmp_out );
-
-while ( $entry = get_tab_entry( $fh_in ) )
-{
-    $record = uclust_tab2biopiece( $entry );
-
-    Maasha::Biopieces::put_record( $record, $out ) if $record->{ 'TYPE' } eq 'Hit';
-}
-
-close $fh_out;
-
-unlink $tmp_out;
-
-Maasha::Biopieces::close_stream( $in );
-Maasha::Biopieces::close_stream( $out );
-
-
-# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-
-sub get_tab_entry
-{
-    # Martin A. Hansen, May 2009.
-
-    # Get the next tabular entry from a filehandle to a USEARCH file.
-
-    my ( $fh,   # filehandle
-       ) = @_;
-
-    # Returns a list
-
-    my ( $line, @fields );
-
-    while ( $line = <$fh> )
-    {
-        next if $line =~ /^#/;
-    
-        @fields = split /\s+/, $line;
-
-        return wantarray ? @fields : \@fields;
-    }
-
-    return undef;
-}
-
-
-sub uclust_tab2biopiece
-{
-    # Martin A. Hansen, May 2009.
-
-    # Get the next USEARCH tabular entry and convert it to
-    # a biopiece record that is returned.
-
-    my ( $entry,    # USEARCH tabular entry
-       ) = @_;
-
-    # Returns a hashref.
-
-    my ( %record );
-
-    $record{ "REC_TYPE" }  = "USEARCH";
-    $record{ "TYPE" }      = $TYPE_HASH->{ $entry->[ 0 ] };
-    $record{ "CLUSTER" }   = $entry->[ 1 ];
-    $record{ "ALIGN_LEN" } = $entry->[ 2 ];
-    $record{ "ID" }        = $entry->[ 3 ];
-    $record{ "STRAND" }    = $entry->[ 4 ];
-    $record{ "Q_BEG" }     = $entry->[ 5 ];
-    $record{ "S_BEG" }     = $entry->[ 6 ];
-    $record{ "CIGAR" }     = $entry->[ 7 ];
-    $record{ "Q_ID" }      = $entry->[ 8 ];
-    $record{ "S_ID" }      = $entry->[ 9 ];
-
-    if ( $record{ 'TYPE' } eq 'Hit' )
-    {
-        $record{ "Q_END" } = $record{ "Q_BEG" } + $record{ "ALIGN_LEN" } - 1;
-        $record{ "S_END" } = $record{ "S_BEG" } + $record{ "ALIGN_LEN" } - 1;
-    }
-
-    return wantarray ? %record : \%record;
-}
-
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
+# Usearch sequences in the stream against a specified database.
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
+require 'maasha/biopieces'
+require 'maasha/fasta'
+require 'maasha/usearch'
 
-BEGIN
-{
-    Maasha::Biopieces::status_set();
-}
+casts = []
+casts << {long: 'program',    short: 'p', type: 'string', mandatory: false, default: 'global', allowed: "local,global", disallowed: nil}
+casts << {long: 'database',   short: 'd', type: 'file!',  mandatory: true,  default: nil,      allowed: nil,            disallowed: nil}
+casts << {long: 'identity',   short: 'i', type: 'float',  mandatory: false, default: nil,      allowed: nil,            disallowed: nil}
+casts << {long: 'e_val',      short: 'e', type: 'float',  mandatory: false, default: nil,      allowed: nil,            disallowed: nil}
+casts << {long: 'cpus',       short: 'c', type: 'uint',   mandatory: false, default: 1,        allowed: nil,            disallowed: "0"}
+casts << {long: 'maxaccepts', short: 'm', type: 'uint',   mandatory: false, default: 0,        allowed: nil,            disallowed: nil}
 
+options = Biopieces.options_parse(ARGV, casts)
 
-END
-{
-    Maasha::Biopieces::status_log();
-}
+if options[:program] == 'global'
+  raise ArgumentError, "--identity be specified for global search" unless options[:identity]
+  raise ArgumentError, "--e_val is invalid for global search" if options[:e_val]
+else
+  raise ArgumentError, "--identity or --e_val must be specified" unless options[:identity] or options[:e_val]
+end
+
+tmpdir  = Biopieces.mktmpdir
+infile  = File.join(tmpdir, "in.fna")
+outfile = File.join(tmpdir, "out.uc")
+
+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
+
+      if record[:SEQ_NAME] and record[:SEQ]
+        fasta_io.puts Seq.new_bp(record).to_fasta
+      end
+    end
+  end
+
+  us = Usearch.new(infile, outfile, options)
+
+  case options[:program]
+  when "local"  then us.usearch_local
+  when "global" then us.usearch_global
+  else raise "no such program #{options[:program]}"
+  end
+
+  us.each_hit do |record|
+    output.puts record
+  end
+end
 
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<