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