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, 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, :match
 113
 114     def initialize(start, stop, strand, match)
 115       @start   = start
 116       @stop    = stop
 117       @strand  = strand
 118       @match   = match
 119     end
 120
 121     def length
 122       @stop - @start + 1
 123     end
 124   end
 125 end
 126
 127 casts = []
 128 casts << {:long=>'pattern',    :short=>'p', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 129 casts << {:long=>'pattern_in', :short=>'P', :type=>'file!',  :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 130 casts << {:long=>'inline',     :short=>'i', :type=>'flag',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 131 casts << {:long=>'overlap',    :short=>'o', :type=>'flag',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 132 casts << {:long=>'comp',       :short=>'c', :type=>'flag',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
 133 casts << {:long=>'max_hits',   :short=>'h', :type=>'uint',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>'0'}
 134 casts << {:long=>'max_misses', :short=>'m', :type=>'uint',   :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>'0'}
 135
 136 options = Biopieces.options_parse(ARGV, casts)
 137
 138 unless options[:pattern] or options[:pattern_in]
 139   raise ArgumentError, "--pattern or --pattern_in must be specified"
 140 end
 141
 142 patterns = []
 143
 144 if options[:pattern_in]
 145   File.open(options[:pattern_in], 'r') do |ios|
 146     ios.each_line { |l| patterns << l.chomp }
 147   end
 148 else
 149   patterns << options[:pattern]
 150 end
 151
 152 raise ArgumentError, "no patterns found" if patterns.empty?
 153
 154 tmp_dir  = Biopieces.mktmpdir
 155 tmp_file = File.join(tmp_dir, "tmp.stream")
 156 in_file  = File.join(tmp_dir, "in.fna")
 157
 158 seq_name_count = 0
 159 seq_name_hash  = {}
 160 seq_type       = nil
 161
 162 Biopieces.open(options[:stream_in], tmp_file) do |input, output|
 163   Fasta.open(in_file, mode="w") do |fasta_io|
 164     input.each_record do |record|
 165       if record[:SEQ_NAME]
 166         seq_name_hash[seq_name_count] = record[:SEQ_NAME]
 167         record[:SEQ_NAME] = seq_name_count
 168         seq_name_count += 1
 169
 170         seq = Seq.new_bp(record)
 171
 172         if seq_type.nil?
 173           seq_type = seq.type_guess
 174         end
 175
 176         fasta_io.puts seq.to_fasta
 177       end
 178
 179       output.puts record
 180     end
 181   end
 182 end
 183
 184 options[:seq_type] = seq_type
 185
 186 patscan = Patscan.new(tmp_dir, in_file, patterns)
 187 patscan.run(options)
 188 results = patscan.results_parse
 189
 190 Biopieces.open(tmp_file, options[:stream_out]) do |input, output|
 191   input.each_record do |record|
 192     if record[:SEQ_NAME]
 193       key = record[:SEQ_NAME].to_i
 194       record[:SEQ_NAME] = seq_name_hash[key]
 195
 196       if options[:inline]
 197         if results[key]
 198           results[key].each do |result|
 199             record[:PATTERN]   = options[:pattern]
 200             record[:MATCH]     = result.match
 201             record[:S_BEG]     = result.start
 202             record[:S_END]     = result.stop
 203             record[:STRAND]    = result.strand
 204             record[:MATCH_LEN] = result.length
 205
 206             output.puts record
 207           end
 208         else
 209           output.puts record
 210         end
 211       else
 212         output.puts record
 213
 214         new_record = {}
 215
 216         results[key].each do |result|
 217           new_record[:REC_TYPE]  = "PATSCAN"
 218           new_record[:S_ID]      = record[:SEQ_NAME]
 219           new_record[:Q_ID]      = options[:pattern]
 220           new_record[:MATCH]     = result.match
 221           new_record[:S_BEG]     = result.start
 222           new_record[:S_END]     = result.stop
 223           new_record[:STRAND]    = result.strand
 224           new_record[:SCORE]     = 100
 225           new_record[:MATCH_LEN] = result.length
 226
 227           output.puts new_record
 228         end
 229       end
 230     else
 231       output.puts record
 232     end
 233   end
 234 end
 235
 236
 237 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 238
 239
 240 __END__