code_ruby/lib/maasha/seq/patscan.rb

   1 # Copyright (C) 2007-2011 Martin A. Hansen.
   2
   3 # This program is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU General Public License
   5 # as published by the Free Software Foundation; either version 2
   6 # of the License, or (at your option) any later version.
   7
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program; if not, write to the Free Software
  15 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16
  17 # http://www.gnu.org/copyleft/gpl.html
  18
  19 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  20
  21 # This software is part of the Biopieces framework (www.biopieces.org).
  22
  23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  24
  25 require 'maasha/filesys'
  26 require 'open3'
  27
  28 # Module containing code to locate nucleotide patterns in sequences allowing for
  29 # ambiguity codes and a given maximum mismatches, insertions, and deletions. The
  30 # pattern match engine is based on scan_for_matches.
  31 module Patscan
  32   # ------------------------------------------------------------------------------
  33   #   str.scan(pattern[, max_mismatches [, max_insertions [,max_deletions]]])
  34   #   -> Array
  35   #   str.scan(pattern[, max_mismatches [, max_insertions [,max_deletions]]]) { |match|
  36   #     block
  37   #   }
  38   #   -> Match
  39   #
  40   # ------------------------------------------------------------------------------
  41   # Method to iterate through a sequence to locate pattern matches starting
  42   # allowing for a maximum number of mismatches, insertions, and deletions.
  43   # Matches found in block context return the Match object. Otherwise matches are
  44   # returned in an Array.
  45   def scan(pattern, max_mismatches = 0, max_insertions = 0, max_deletions = 0)
  46     @pattern = pattern + "[#{max_mismatches},#{max_insertions},#{max_deletions}]"
  47     matches  = []
  48
  49     pattern_file = Filesys.tmpfile
  50
  51     File.open(pattern_file, mode = 'w') do |ios|
  52       ios.puts @pattern
  53     end
  54
  55     cmd = "scan_for_matches #{pattern_file}"
  56
  57     Open3.popen3(cmd) do |stdin, stdout, stderr|
  58       #raise SeqError, "scan_for_matches failed: #{stderr.gets}" unless stderr.empty?
  59
  60       stdin << ">dummy\n#{self.seq}\n"
  61
  62       stdin.close
  63
  64       while block = stdout.gets($/ + ">")
  65         if block =~ />dummy:\[(\d+),(\d+)\]\n(.+)/
  66           start  = $1.to_i
  67           stop   = $2.to_i
  68           match  = $3.delete(" ")
  69           length = stop - start
  70
  71           if block_given?
  72             yield Match.new(start, length, match)
  73           else
  74             matches << Match.new(start, length, match)
  75           end
  76         end
  77       end
  78     end
  79
  80     matches unless block_given?
  81   end
  82
  83   class Match
  84     attr_reader :pos, :length, :match
  85
  86     def initialize(pos, length, match)
  87       @pos    = pos
  88       @length = length
  89       @match  = match
  90     end
  91   end
  92 end
  93
  94 __END__
  95
  96
  97 # Eeeeiiiii - this one is slower than scan_for_matches!
  98
  99
 100 # Module containing code to locate nucleotide patterns in sequences allowing for
 101 # ambiguity codes and a given maximum edit distance. Pattern match engine is based
 102 # on nrgrep.
 103 module Patscan
 104   # ------------------------------------------------------------------------------
 105   #   str.scan(pattern[, max_edit_distance])
 106   #   -> Array
 107   #   str.scan(pattern[, max_edit_distance]) { |match|
 108   #     block
 109   #   }
 110   #   -> Match
 111   #
 112   # ------------------------------------------------------------------------------
 113   # Method to iterate through a sequence to locate pattern matches starting
 114   # from a given position and allowing for a maximum edit distance.
 115   # Matches found in block context return the Match object. Otherwise matches are
 116   # returned in an Array.
 117   def scan(pattern, edit_distance = 0)
 118     @pattern = pattern.upcase
 119     matches  = []
 120
 121     pattern_disambiguate
 122
 123     cmd = "nrgrep_coords"
 124     cmd << " -i"
 125     cmd << " -k #{edit_distance}s" if edit_distance > 0
 126     cmd << " #{@pattern}"
 127
 128     Open3.popen3(cmd) do |stdin, stdout, stderr|
 129       stdin << self.seq
 130
 131       stdin.close
 132
 133       stdout.each_line do |line|
 134         if line =~ /\[(\d+), (\d+)\]: (.+)/
 135           start = $1.to_i
 136           stop  = $2.to_i
 137           match = $3
 138           length = stop - start
 139
 140           if block_given?
 141             yield Match.new(start, length, match)
 142           else
 143             matches << Match.new(start, length, match)
 144           end
 145         end
 146       end
 147     end
 148
 149     matches unless block_given?
 150   end
 151
 152   private
 153
 154   def pattern_disambiguate
 155     case self.type
 156     when :protein
 157       @pattern.gsub!('J', '[IFVLWMAGCY]')
 158       @pattern.gsub!('O', '[TSHEDQNKR]')
 159       @pattern.gsub!('B', '[DN]')
 160       @pattern.gsub!('Z', '[EQ]')
 161       @pattern.gsub!('X', '.')
 162     when :dna
 163       @pattern.gsub!('R', '[AG]')
 164       @pattern.gsub!('Y', '[CT]')
 165       @pattern.gsub!('S', '[GC]')
 166       @pattern.gsub!('W', '[AT]')
 167       @pattern.gsub!('M', '[AC]')
 168       @pattern.gsub!('K', '[GT]')
 169       @pattern.gsub!('V', '[ACG]')
 170       @pattern.gsub!('H', '[ACT]')
 171       @pattern.gsub!('D', '[AGT]')
 172       @pattern.gsub!('B', '[CGT]')
 173       @pattern.gsub!('N', '.')
 174     when :rna
 175     else
 176       raise SeqError "unknown sequence type: #{self.type}"
 177     end
 178   end
 179
 180   class Match
 181     attr_reader :pos, :length, :match
 182
 183     def initialize(pos, length, match)
 184       @pos    = pos
 185       @length = length
 186       @match  = match
 187     end
 188   end
 189 end
 190