bp_bin/kmer_stats

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2012 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 # Replace values to a specied key for each record in the stream.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31 require 'maasha/biopieces'
  32 require 'maasha/seq'
  33 require 'inline'
  34
  35 # Class containing methods to locate over-represented kmers in sequences.
  36 class KmerStats
  37   NUC_ALPH_SIZE = 4
  38   BYTES_IN_INT  = 4
  39
  40   # Method to initialize a KmerStats object.
  41   def initialize(kmer_size)
  42     @kmer_size = kmer_size
  43
  44     raise ArgumentError, "kmer_size: #{kmer_size} > 16" if kmer_size > 16
  45
  46     @kmer_tot        = 0
  47     @freq_tot        = 0
  48     @kmer_index_size = NUC_ALPH_SIZE ** kmer_size
  49     @freq_index_size = NUC_ALPH_SIZE
  50     @kmer_index      = "\0" * @kmer_index_size * BYTES_IN_INT
  51     @freq_index      = "\0" * @freq_index_size * BYTES_IN_INT
  52   end
  53
  54   # Method to add a sequence to the KmerStats.
  55   def add(seq)
  56     @kmer_tot += add_C(seq, seq.length, @kmer_size, @kmer_index, @freq_index)
  57     @freq_tot += seq.length
  58   end
  59
  60   # Method for iterating the KmerStats and return.
  61   # each kmer as a binary value.
  62   def each
  63     (0 ... @kmer_index_size).each do |bin|
  64       yield bin
  65     end
  66   end
  67
  68   # Method to determine the probability of an observed binary kmer.
  69   def p_kmer(bin)
  70     p_kmer_C(bin, @kmer_tot, @kmer_index)
  71   end
  72
  73   # Method to determine the probablity of a theoretical binary kmer.
  74   def p_freq(bin)
  75     p_freq_C(bin, @freq_tot, @kmer_size, @freq_index)
  76   end
  77
  78   # Method that converts a binary kmer to DNA.
  79   def kmer(bin)
  80     dna  = ""
  81     mask = 3
  82
  83     (0 ... @kmer_size).each do
  84       case bin & mask
  85       when 1 then dna << "T"
  86       when 2 then dna << "C"
  87       when 3 then dna << "G"
  88       else
  89         dna << "A"
  90       end
  91
  92       bin >>= 2
  93     end
  94
  95     dna.reverse
  96   end
  97
  98   # Method that returns the count of a given binary kmer.
  99   def count(bin)
 100     count_C(bin, @kmer_index)
 101   end
 102
 103   private
 104
 105   inline do |builder|
 106     builder.prefix %{
 107       unsigned char nuc2bin[256] = {
 108         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 109         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 110         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 111         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 112         0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0,
 113         0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 114         0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0,
 115         0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 116         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 117         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 118         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 119         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 120         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 121         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 122         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 123         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 124       };
 125     }
 126
 127     builder.c %{
 128       VALUE add_C(
 129         VALUE _seq,
 130         VALUE _seq_len,
 131         VALUE _kmer_size,
 132         VALUE _kmer_index,
 133         VALUE _freq_index
 134       )
 135       {
 136         unsigned char *seq        = (unsigned char *) StringValuePtr(_seq);
 137         unsigned int   seq_len    = NUM2UINT(_seq_len);
 138         unsigned int   kmer_size  = NUM2UINT(_kmer_size);
 139         unsigned int  *kmer_index = (unsigned int *) StringValuePtr(_kmer_index);
 140         unsigned int  *freq_index = (unsigned int *) StringValuePtr(_freq_index);
 141
 142         unsigned int mask = (1 << (2 * kmer_size)) - 1;
 143         unsigned int bin  = 0;
 144         unsigned int val  = 0;
 145         unsigned int i    = 0;
 146
 147         for (i = 0; i < kmer_size; i++)
 148         {
 149           bin <<= 2;
 150           val = nuc2bin[seq[i]];
 151           bin |= val;
 152           freq_index[val]++;
 153         }
 154
 155         kmer_index[bin]++;
 156
 157         while (i < seq_len)
 158         {
 159           bin <<= 2;
 160           val = nuc2bin[seq[i]];
 161           bin |= val;
 162
 163           kmer_index[bin & mask]++;
 164           freq_index[val]++;
 165
 166           i++;
 167         }
 168
 169         return UINT2NUM(i - kmer_size + 1);
 170       }
 171     }
 172
 173     builder.c %{
 174       VALUE freq_count_C(
 175         VALUE _seq,
 176         VALUE _seq_len,
 177         VALUE _index
 178       )
 179       {
 180         unsigned char *seq     = (unsigned char *) StringValuePtr(_seq);
 181         unsigned int   seq_len = NUM2UINT(_seq_len);
 182         unsigned int  *index   = (unsigned int *) StringValuePtr(_index);
 183
 184         unsigned int i = 0;
 185
 186         for (i = 0; i < seq_len; i++)
 187         {
 188           index[seq[i]]++;
 189         }
 190
 191         return UINT2NUM(i);
 192       }
 193     }
 194
 195     builder.c %{
 196       VALUE p_kmer_C(
 197         VALUE _bin,
 198         VALUE _kmer_tot,
 199         VALUE _kmer_index
 200       )
 201       {
 202         unsigned int  bin        = NUM2UINT(_bin);
 203         unsigned int  kmer_tot   = NUM2UINT(_kmer_tot);
 204         unsigned int *kmer_index = (unsigned int *) StringValuePtr(_kmer_index);
 205
 206         double p_kmer = 0.0;
 207
 208         p_kmer = (double) kmer_index[bin] / kmer_tot;
 209
 210         return DBL2NUM(p_kmer);
 211       }
 212     }
 213
 214     builder.c %{
 215       VALUE p_freq_C(
 216         VALUE _bin,
 217         VALUE _freq_tot,
 218         VALUE _kmer_size,
 219         VALUE _freq_index
 220       )
 221       {
 222         unsigned int  bin        = NUM2UINT(_bin);
 223         unsigned int  freq_tot   = NUM2UINT(_freq_tot);
 224         unsigned int  kmer_size  = NUM2UINT(_kmer_size);
 225         unsigned int *freq_index = (unsigned int *) StringValuePtr(_freq_index);
 226
 227         unsigned int mask   = 3;
 228         unsigned int i      = 0;
 229         unsigned int c      = 0;
 230         double       p_freq = 1.0;
 231
 232         for (i = 0; i < kmer_size; i++)
 233         {
 234           p_freq *= (double) freq_index[bin & mask] / freq_tot;
 235
 236           bin >>= 2;
 237         }
 238
 239         return DBL2NUM(p_freq);
 240       }
 241     }
 242
 243     builder.c %{
 244       VALUE count_C(
 245         VALUE _bin,
 246         VALUE _kmer_index
 247       )
 248       {
 249         unsigned int  bin        = NUM2UINT(_bin);
 250         unsigned int *kmer_index = (unsigned int *) StringValuePtr(_kmer_index);
 251
 252         unsigned int count = kmer_index[bin];
 253
 254         return UINT2NUM(count);
 255       }
 256     }
 257   end
 258 end
 259
 260 casts = []
 261 casts << {:long=>'size',      :short=>'s', :type=>'uint',  :mandatory=>true,  :default=>10, :allowed=>nil, :disallowed=>"0"}
 262 casts << {:long=>'min_ratio', :short=>'m', :type=>'float', :mandatory=>false, :default=>10, :allowed=>nil, :disallowed=>nil}
 263
 264 options = Biopieces.options_parse(ARGV, casts)
 265
 266 kmer_stats = KmerStats.new(options[:size])
 267
 268 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
 269   input.each_record do |record|
 270     output.puts record
 271
 272     if record[:SEQ]
 273       kmer_stats.add(record[:SEQ])
 274     end
 275   end
 276
 277   kmer_stats.each do |bin|
 278     p_kmer = kmer_stats.p_kmer(bin)
 279
 280     next if p_kmer == 0.0
 281
 282     p_freq = kmer_stats.p_freq(bin)
 283
 284     ratio = (p_kmer / p_freq).round(2)
 285
 286     if ratio >= options[:min_ratio]
 287       kmer  = kmer_stats.kmer(bin)
 288       count = kmer_stats.count(bin)
 289
 290       record = {
 291         :REC_TYPE => "KMER_STATS",
 292         :P_KMER   => p_kmer,
 293         :KMER     => kmer,
 294         :COUNT    => count,
 295         :RATIO    => ratio
 296       }
 297
 298       output.puts record
 299     end
 300   end
 301 end
 302
 303 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 304
 305
 306 __END__