bp_bin/bwa_seq

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2012 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 # Use BWA to map sequences in the stream against a specified index.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31 require 'maasha/biopieces'
  32 require 'maasha/seq'
  33 require 'maasha/fasta'
  34 require 'maasha/fastq'
  35 require 'maasha/sam'
  36
  37 casts = []
  38 casts << {:long=>'index_name', :short=>'i', :type=>'string', :mandatory=>true,  :default=>nil, :allowed=>nil,       :disallowed=>nil}
  39 casts << {:long=>'mismatches', :short=>'m', :type=>'uint',   :mandatory=>false, :default=>0,   :allowed=>"0,1,2,3", :disallowed=>nil}
  40 casts << {:long=>'max_hits',   :short=>'h', :type=>'uint',   :mandatory=>false, :default=>30,  :allowed=>nil,       :disallowed=>'0'}
  41 casts << {:long=>'cpus',       :short=>'c', :type=>'uint',   :mandatory=>false, :default=>1,   :allowed=>nil,       :disallowed=>'0'}
  42
  43 options = Biopieces.options_parse(ARGV, casts)
  44
  45 tmp_dir = Biopieces.mktmpdir
  46 tmp_fq  = File.join(tmp_dir, "bwa.fq")
  47 tmp_sai = File.join(tmp_dir, "bwa.sai")
  48 tmp_sam = File.join(tmp_dir, "bwa.sam")
  49
  50 # Class containing methods for executing BWA.
  51 class BWA
  52   def initialize(fq_file, sai_file, sam_file, options)
  53     @fq_file  = fq_file
  54     @sai_file = sai_file
  55     @sam_file = sam_file
  56     @options  = options
  57     @command  = []
  58   end
  59
  60   # Method to run bwa aln.
  61   def aln
  62     @command << "bwa aln"
  63     @command << "-t #{@options[:cpus]}"
  64     @command << "-R #{@options[:max_hits]}"
  65     @command << "-f #{@sai_file}"
  66     @command << @options[:index_name]
  67     @command << @fq_file
  68
  69     execute
  70   end
  71
  72   # Method to run bwa samse.
  73   def samse
  74     @command << "bwa samse"
  75     @command << "-f #{@sam_file}"
  76     @command << @options[:index_name]
  77     @command << @sai_file
  78     @command << @fq_file
  79
  80     execute
  81   end
  82
  83   private
  84
  85   # Method to execute a command using a system() call.
  86   # The command is composed of bits from the @command variable.
  87   def execute
  88     @command.unshift "nice -n 19"
  89     @command << "> /dev/null 2>&1" unless @options[:verbose]
  90
  91     command = @command.join(" ")
  92     $stderr.puts "Running command: #{command}" if @options[:verbose]
  93     system(command)
  94     raise "Command failed: #{command}" unless $?.success?
  95
  96     @command = []
  97   end
  98 end
  99
 100 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
 101   Fastq.open(tmp_fq, 'w') do |io_fq|
 102     input.each_record do |record|
 103       output.puts record
 104
 105       if record[:SEQ_NAME] and record[:SEQ] and record[:SCORES]
 106         entry = Seq.new_bp(record)
 107
 108         io_fq.puts entry.to_fastq
 109       end
 110     end
 111   end
 112
 113   bwa = BWA.new(tmp_fq, tmp_sai, tmp_sam, options)
 114   bwa.aln
 115   bwa.samse
 116
 117   Sam.open(tmp_sam, 'r') do |io_sam|
 118     io_sam.each do |entry|
 119       output.puts Sam.to_bp(entry) unless entry[:RNAME] == '*'
 120     end
 121   end
 122 end
 123