bp_bin/patscan_seq

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2011 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 Patscan on sequences in the stream.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31
  32 require 'maasha/biopieces'
  33 require 'maasha/fasta'
  34 require 'maasha/seq'
  35
  36 # Class with methods to execute Patscan (a.k.a. scan_for_matches).
  37 class Patscan
  38   # Method to initialize a Patscan object.
  39   def initialize(pat_file, in_file, out_file, pattern)
  40     @pat_file = pat_file
  41     @in_file  = in_file
  42     @out_file = out_file
  43     @pattern  = pattern
  44
  45     pattern_save
  46   end
  47
  48   # Method to run Patscan
  49   def run(comp = nil, type = nil, max_hit = nil, max_mis = nil, verbose = nil)
  50     args = []
  51     args << "scan_for_matches"
  52     args << "-c"            if comp
  53     args << "-p"            if type == :protein
  54     args << "-n #{max_mis}" if max_mis
  55     args << "-m #{max_hit}" if max_hit
  56     args << @pat_file
  57     args << "< #{@in_file}"
  58     args << "> #{@out_file}"
  59                 command = args.join(" ")
  60     $stderr.puts command if verbose
  61     system(command)
  62     raise "Command failed: #{command}" unless $?.success?
  63   end
  64
  65   # Method to parse Patscan results and return
  66   # these in a hash with ID as key and a list
  67   # of hits as value.
  68   def results_parse
  69     results = Hash.new { |h, k| h[k] = [] }
  70
  71     Fasta.open(@out_file, 'r') do |ios|
  72       ios.each do |entry|
  73         if entry.seq_name =~ /([^:]+):\[(\d+),(\d+)\]/
  74           id    = $1.to_i
  75           start = $2.to_i - 1
  76           stop  = $3.to_i - 1
  77
  78           if stop > start
  79             strand = '+'
  80           else
  81             start, stop = stop, start
  82             strand = '-'
  83           end
  84
  85           results[id.to_i] << Match.new(start, stop, strand, entry.seq)
  86         else
  87           raise "Failed to parse seq_name: #{entry.seq_name} in patscan result"
  88         end
  89       end
  90     end
  91
  92     results
  93   end
  94
  95   private
  96
  97   # Method to save a patscan pattern to a file.
  98   def pattern_save
  99     File.open(@pat_file, "w") do |ios|
 100       ios.puts @pattern
 101     end
 102   end
 103
 104   # Subclass to define Patscan hits.
 105   class Match
 106     attr_reader :start, :stop, :strand, :match
 107
 108     def initialize(start, stop, strand, match)
 109       @start   = start
 110       @stop    = stop
 111       @strand  = strand
 112       @match   = match
 113     end
 114
 115     def length
 116       @stop - @start + 1
 117     end
 118   end
 119 end
 120
 121 casts = []
 122 casts << {:long=>'pattern',    :short=>'p', :type=>'string', :mandatory=>true,  :default=>nil, :allowed=>nil, :disallowed=>nil}
 123 casts << {:long=>'inline',     :short=>'i', :type=>'flag',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 124 casts << {:long=>'comp',       :short=>'c', :type=>'flag',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 125 casts << {:long=>'max_hits',   :short=>'h', :type=>'uint',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>'0'}
 126 casts << {:long=>'max_misses', :short=>'m', :type=>'uint',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>'0'}
 127
 128 options = Biopieces.options_parse(ARGV, casts)
 129
 130 tmp_dir  = Biopieces.mktmpdir
 131 tmp_file = File.join(tmp_dir, "tmp.stream")
 132 pat_file = File.join(tmp_dir, "pat.txt")
 133 in_file  = File.join(tmp_dir, "in.fna")
 134 out_file = File.join(tmp_dir, "out.fna")
 135 out_file = File.join(tmp_dir, "out.fna")
 136
 137 seq_name_count = 0
 138 seq_name_hash  = {}
 139 seq_type       = nil
 140
 141 Biopieces.open(options[:stream_in], tmp_file) do |input, output|
 142   Fasta.open(in_file, mode="w") do |fasta_io|
 143     input.each_record do |record|
 144       if record[:SEQ_NAME]
 145         seq_name_hash[seq_name_count] = record[:SEQ_NAME]
 146         record[:SEQ_NAME] = seq_name_count
 147         seq_name_count += 1
 148
 149         seq = Seq.new_bp(record)
 150
 151         if seq_type.nil?
 152           seq_type = seq.type_guess
 153         end
 154
 155         fasta_io.puts seq.to_fasta
 156       end
 157
 158       output.puts record
 159     end
 160   end
 161 end
 162
 163 patscan = Patscan.new(pat_file, in_file, out_file, options[:pattern])
 164 patscan.run(options[:comp], seq_type, options[:max_hits], options[:max_misses], options[:verbose])
 165 results = patscan.results_parse
 166
 167 Biopieces.open(tmp_file, options[:stream_out]) do |input, output|
 168   input.each_record do |record|
 169     if record[:SEQ_NAME]
 170       key = record[:SEQ_NAME].to_i
 171       record[:SEQ_NAME] = seq_name_hash[key]
 172
 173       if options[:inline]
 174         if results[key]
 175           results[key].each do |result|
 176             record[:PATTERN]   = options[:pattern]
 177             record[:MATCH]     = result.match
 178             record[:S_BEG]     = result.start
 179             record[:S_END]     = result.stop
 180             record[:STRAND]    = result.strand
 181             record[:MATCH_LEN] = result.length
 182
 183             output.puts record
 184           end
 185         else
 186           output.puts record
 187         end
 188       else
 189         output.puts record
 190
 191         new_record = {}
 192
 193         results[key].each do |result|
 194           new_record[:REC_TYPE]  = "PATSCAN"
 195           new_record[:S_ID]      = record[:SEQ_NAME]
 196           new_record[:Q_ID]      = options[:pattern]
 197           new_record[:MATCH]     = result.match
 198           new_record[:S_BEG]     = result.start
 199           new_record[:S_END]     = result.stop
 200           new_record[:STRAND]    = result.strand
 201           new_record[:SCORE]     = 100
 202           new_record[:MATCH_LEN] = result.length
 203
 204           output.puts new_record
 205         end
 206       end
 207     else
 208       output.puts record
 209     end
 210   end
 211 end
 212
 213
 214 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 215
 216
 217 __END__