X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bp_bin%2Fbwa_seq;h=fd807821dac3dfe3728a8b34dae9f7ee1ec25fc3;hb=5de6112b70b59420b245ce636a8b2e3c90acbe00;hp=c712b8f86a17c2de2bebdc515e13245e5d9b8056;hpb=b4039cd22036599a2ea32211ff50f20cc9345eab;p=biopieces.git

diff --git a/bp_bin/bwa_seq b/bp_bin/bwa_seq
index c712b8f..fd80782 100755
--- a/bp_bin/bwa_seq
+++ b/bp_bin/bwa_seq
@@ -1,6 +1,6 @@
-#!/usr/bin/env perl
+#!/usr/bin/env ruby
 
-# Copyright (C) 2007-2009 Martin A. Hansen.
+# Copyright (C) 2007-2012 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,121 +18,106 @@
 
 # http://www.gnu.org/copyleft/gpl.html
 
-
-# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-# Use BWA to map sequences in the stream against a specified genome or index.
-
-# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-
-use warnings;
-use strict;
-use Data::Dumper;
-use Maasha::Biopieces;
-use Maasha::Common;
-use Maasha::Fastq;
-use Maasha::SAM;
-
-
-# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-
-my ( $options, $in, $out, $index, $tmp_dir, $tmp_fq, $tmp_sai, $tmp_sam, $fh_in, $fh_out, $record, $entry, $line, @fields, $type, @args, $arg );
-
-$options = Maasha::Biopieces::parse_options(
-    [
-        { long => 'genome',      short => 'g', type => 'genome', mandatory => 'no', default => undef, allowed => undef,     disallowed => undef },
-        { long => 'index_name',  short => 'i', type => 'string', mandatory => 'no', default => undef, allowed => undef,     disallowed => undef },
-        { long => 'mismatches',  short => 'm', type => 'uint',   mandatory => 'no', default => 0,     allowed => "0,1,2,3", disallowed => undef },
-        { long => 'max_hits',    short => 'h', type => 'uint',   mandatory => 'no', default => undef, allowed => undef,     disallowed => 0     },
-        { long => 'cpus',        short => 'c', type => 'uint',   mandatory => 'no', default => 1,     allowed => undef,     disallowed => 0     },
-    ]   
-);
-
-Maasha::Common::error( qq(both --index_name and --genome specified) ) if     $options->{ "genome" } and     $options->{ "index_name" };
-Maasha::Common::error( qq(no --index_name or --genome specified) )    if not $options->{ "genome" } and not $options->{ "index_name" };
-
-$in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
-$out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
-
-if ( defined $options->{ 'genome' } ) {
-    $index = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/bowtie/$options->{ 'genome' }";
-} elsif (defined $options->{ 'index_name' } ) {
-    $index = $options->{ 'index_name' };
-}
-
-$tmp_dir = Maasha::Biopieces::get_tmpdir();
-$tmp_fq  = "$tmp_dir/bwa.fq";
-$tmp_sai = "$tmp_dir/bwa.fq.sai";
-$tmp_sam = "$tmp_dir/bwa.fq.sam";
-
-$fh_out = Maasha::Filesys::file_write_open( $tmp_fq );
-
-while ( $record = Maasha::Biopieces::get_record( $in ) ) 
-{
-    if ( $entry = Maasha::Fastq::biopiece2fastq( $record ) ) {
-        Maasha::Fastq::put_entry( $entry, $fh_out );
-    }
-
-    Maasha::Biopieces::put_record( $record, $out );
-}
-
-close $fh_out;
-
-push @args, "-t $options->{ 'cpus' }";
-push @args, "-R $options->{ 'max_hits' }" if $options->{ 'max_hits' };
-
-$arg = join " ", @args;
-
-if ( $options->{ 'verbose' } )
-{
-    print STDERR qq(Running: bwa aln $arg $index $tmp_fq > $tmp_sai\n);
-    Maasha::Common::run( "bwa", "aln $arg $index $tmp_fq > $tmp_sai" );
-    print STDERR qq(Running: bwa samse $index $tmp_sai $tmp_fq > $tmp_sam\n);
-    Maasha::Common::run( "bwa", "samse $index $tmp_sai $tmp_fq > $tmp_sam" )
-}
-else
-{
-    Maasha::Common::run( "bwa", "aln $arg $index $tmp_fq       > $tmp_sai 2> /dev/null" );
-    Maasha::Common::run( "bwa", "samse $index $tmp_sai $tmp_fq > $tmp_sam 2> /dev/null" );
-}
-
-$fh_in = Maasha::Filesys::file_read_open( $tmp_sam );
-
-while ( $entry = Maasha::SAM::get_entry( $fh_in ) )
-{
-    if ( $record = Maasha::SAM::sam2biopiece( $entry ) ) {
-        Maasha::Biopieces::put_record( $record, $out );
-    }
-}
-
-close $fh_in;
-
-unlink $tmp_fq;
-unlink $tmp_sai;
-unlink $tmp_sam;
-
-Maasha::Biopieces::close_stream( $in );
-Maasha::Biopieces::close_stream( $out );
-
-
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
+# This program is part of the Biopieces framework (www.biopieces.org).
 
-BEGIN
-{
-    Maasha::Biopieces::status_set();
-}
-
-
-END
-{
-    Maasha::Biopieces::status_log();
-}
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
+# Use BWA to map sequences in the stream against a specified index.
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
+require 'maasha/biopieces'
+require 'maasha/seq'
+require 'maasha/fasta'
+require 'maasha/fastq'
+require 'maasha/sam'
+
+casts = []
+casts << {:long=>'index_name', :short=>'i', :type=>'string', :mandatory=>true,  :default=>nil, :allowed=>nil,       :disallowed=>nil}
+casts << {:long=>'mismatches', :short=>'m', :type=>'uint',   :mandatory=>false, :default=>0,   :allowed=>"0,1,2,3", :disallowed=>nil}
+casts << {:long=>'max_hits',   :short=>'h', :type=>'uint',   :mandatory=>false, :default=>30,  :allowed=>nil,       :disallowed=>'0'}
+casts << {:long=>'cpus',       :short=>'c', :type=>'uint',   :mandatory=>false, :default=>1,   :allowed=>nil,       :disallowed=>'0'}
+
+options = Biopieces.options_parse(ARGV, casts)
+
+tmp_dir = Biopieces.mktmpdir
+tmp_fq  = File.join(tmp_dir, "bwa.fq")
+tmp_sai = File.join(tmp_dir, "bwa.sai")
+tmp_sam = File.join(tmp_dir, "bwa.sam")
+
+# Class containing methods for executing BWA.
+class BWA
+  def initialize(fq_file, sai_file, sam_file, options)
+    @fq_file  = fq_file
+    @sai_file = sai_file
+    @sam_file = sam_file
+    @options  = options
+    @command  = []
+  end
+
+  # Method to run bwa aln.
+  def aln
+    @command << "bwa aln"
+    @command << "-t #{@options[:cpus]}"
+    @command << "-R #{@options[:max_hits]}"
+    @command << "-f #{@sai_file}"
+    @command << @options[:index_name]
+    @command << @fq_file
+
+    execute
+  end
+
+  # Method to run bwa samse.
+  def samse
+    @command << "bwa samse"
+    @command << "-f #{@sam_file}"
+    @command << @options[:index_name]
+    @command << @sai_file
+    @command << @fq_file
+
+    execute
+  end
+
+  private
+
+  # Method to execute a command using a system() call.
+  # The command is composed of bits from the @command variable.
+  def execute
+    @command.unshift "nice -n 19"
+    @command << "> /dev/null 2>&1" unless @options[:verbose]
+
+    command = @command.join(" ")
+    $stderr.puts "Running command: #{command}" if @options[:verbose]
+    system(command)
+    raise "Command failed: #{command}" unless $?.success?
+
+    @command = []
+  end
+end
+
+Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
+  Fastq.open(tmp_fq, 'w') do |io_fq|
+    input.each_record do |record|
+      output.puts record
+
+      if record[:SEQ_NAME] and record[:SEQ] and record[:SCORES]
+        entry = Seq.new_bp(record)
+
+        io_fq.puts entry.to_fastq
+      end
+    end
+  end
+
+  bwa = BWA.new(tmp_fq, tmp_sai, tmp_sam, options)
+  bwa.aln
+  bwa.samse
+
+  Sam.open(tmp_sam, 'r') do |io_sam|
+    io_sam.each do |entry|
+      output.puts Sam.to_bp(entry) unless entry[:RNAME] == '*'
+    end
+  end
+end
 
-__END__