code_ruby/lib/maasha/seq/ambiguity.rb

   1 # Copyright (C) 2013 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 'inline'
  26 #autoload :NArray, 'narray'
  27
  28 module Ambiguity
  29   # IUPAC alphabet and binary encoding of the same
  30   # http://en.wikipedia.org/wiki/Nucleic_acid_notation
  31
  32   AMBIGUITY_STR = "ACGTUWSMKRYBDHVNacgtuwsmkrybdhvn"
  33   AMBIGUITY_BIN = "\x08\x04\x02\x01\x01\x09\x06\x0c\x03\x0a\x05\x07\x0b\x0d\x0e\x0f\x08\x04\x02\x01\x01\x09\x06\x0c\x03\x0a\x05\x07\x0b\x0d\x0e\x0f"
  34
  35   # Class method to convert a sequence string to a bit string
  36   # where the bit positions in each char corresponds to the following:
  37   # A = 1000
  38   # C = 0100
  39   # G = 0010
  40   # T = 0001
  41   # And ambiguity codes are expressed using similar bit fields.
  42   def self.to_bin(seq)
  43     seq.tr(AMBIGUITY_STR, AMBIGUITY_BIN)
  44   end
  45
  46   # Class method to convert a bit string to a NArray.
  47   def self.to_na(seq)
  48     NArray.to_na(self.to_bin(seq), 'byte')
  49   end
  50
  51   # Class method to calculate the Hamming Distance between
  52   # two bit fields encoding in NArrays.
  53   def self.hamming_distance(seq1, seq2)
  54     (self.to_na(seq1) & self.to_na(seq2)).eq(0).sum
  55   end
  56
  57   def add_ambiguity_macro inline_builder
  58     # Macro for matching nucleotides including ambiguity codes.
  59     inline_builder.prefix %{
  60       #define MATCH(A,B) ((bitmap[A] & bitmap[B]) != 0)
  61     }
  62
  63     # Bitmap for matching nucleotides including ambiguity codes.
  64     # For each value bits are set from the left: bit pos 1 for A,
  65     # bit pos 2 for T, bit pos 3 for C, and bit pos 4 for G.
  66     inline_builder.prefix %{
  67       char bitmap[256] = {
  68           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  69           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  70           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  71           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  72           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
  73           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
  74           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
  75           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
  76           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  77           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  78           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  79           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  80           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  81           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  82           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  83           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
  84       };
  85     }
  86   end
  87 end
  88