From a1dea64932bbda4eb2a90de9d3e5d353ebeed295 Mon Sep 17 00:00:00 2001 From: martinahansen Date: Mon, 28 Jun 2010 08:18:36 +0000 Subject: [PATCH] added usearch_seq biopiece git-svn-id: http://biopieces.googlecode.com/svn/trunk@998 74ccb610-7750-0410-82ae-013aeee3265d --- bp_bin/usearch_seq | 225 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 225 insertions(+) create mode 100755 bp_bin/usearch_seq diff --git a/bp_bin/usearch_seq b/bp_bin/usearch_seq new file mode 100755 index 0000000..e563893 --- /dev/null +++ b/bp_bin/usearch_seq @@ -0,0 +1,225 @@ +#!/usr/bin/env perl + +# Copyright (C) 2007-2010 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 +# as published by the Free Software Foundation; either version 2 +# of the License, or (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +# 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; + + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +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; +} + + + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +BEGIN +{ + Maasha::Biopieces::status_set(); +} + + +END +{ + Maasha::Biopieces::status_log(); +} + + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +__END__ -- 2.39.2