code_ruby/lib/maasha/seq/trim.rb

   1 # Copyright (C) 2007-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 require 'maasha/seq/backtrack'
  27
  28 include BackTrack
  29
  30 # Error class for all exceptions to do with Trim.
  31 class TrimError < StandardError; end
  32
  33 # Module containing methods for end trimming sequences with suboptimal quality
  34 # scores.
  35 module Trim
  36   # Method to progressively trim a Seq object sequence from the right end until
  37   # a run of min_len residues with quality scores above min_qual is encountered.
  38   def quality_trim_right(min_qual, min_len = 1)
  39     check_trim_args(min_qual)
  40
  41     pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  42
  43     self.subseq(0, pos)
  44   end
  45
  46   # Method to progressively trim a Seq object sequence from the right end until
  47   # a run of min_len residues with quality scores above min_qual is encountered.
  48   def quality_trim_right!(min_qual, min_len = 1)
  49     check_trim_args(min_qual)
  50
  51     pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  52
  53     self.subseq!(0, pos)
  54   end
  55
  56   # Method to progressively trim a Seq object sequence from the left end until
  57   # a run of min_len residues with quality scores above min_qual is encountered.
  58   def quality_trim_left(min_qual, min_len = 1)
  59     check_trim_args(min_qual)
  60
  61     pos = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  62
  63     self.subseq(pos)
  64   end
  65
  66   # Method to progressively trim a Seq object sequence from the left end until
  67   # a run of min_len residues with quality scores above min_qual is encountered.
  68   def quality_trim_left!(min_qual, min_len = 1)
  69     check_trim_args(min_qual)
  70
  71     pos = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  72
  73     self.subseq!(pos)
  74   end
  75
  76   # Method to progressively trim a Seq object sequence from both ends until a
  77   # run of min_len residues with quality scores above min_qual is encountered.
  78   def quality_trim(min_qual, min_len = 1)
  79     check_trim_args(min_qual)
  80
  81     pos_right = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  82     pos_left  = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  83
  84     pos_left = pos_right if pos_left > pos_right
  85
  86     self.subseq(pos_left, pos_right - pos_left)
  87   end
  88
  89   # Method to progressively trim a Seq object sequence from both ends until a
  90   # run of min_len residues with quality scores above min_qual is encountered.
  91   def quality_trim!(min_qual, min_len = 1)
  92     check_trim_args(min_qual)
  93
  94     pos_right = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  95     pos_left  = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  96
  97     pos_left = pos_right if pos_left > pos_right
  98
  99     self.subseq!(pos_left, pos_right - pos_left)
 100   end
 101
 102   # Method to locate a pattern in a sequence and trim all sequence to the left
 103   # including the pattern.
 104   def patmatch_trim_left!(pattern, options = {})
 105     match = self.patmatch(pattern, options)
 106
 107     if match
 108       stop = self.length
 109
 110       self.seq  = self.seq[match.pos + match.length .. stop]
 111       self.qual = self.qual[match.pos + match.length .. stop] if self.qual
 112
 113       return self
 114     end
 115
 116     nil
 117   end
 118
 119   # Method to locate a pattern in a sequence and trim all sequence to the right
 120   # including the pattern.
 121   def patmatch_trim_right!(pattern, options = {})
 122     match = self.patmatch(pattern, options)
 123
 124     if match
 125       self.seq  = self.seq[0 ... match.pos]
 126       self.qual = self.qual[0 ... match.pos] if self.qual
 127
 128       return self
 129     end
 130
 131     nil
 132   end
 133
 134   private
 135
 136   # Method to check the arguments for trimming and raise on bad sequence, qualities,
 137   # and min_qual.
 138   def check_trim_args(min_qual)
 139     raise TrimError, "no sequence"      if self.seq.nil?
 140     raise TrimError, "no quality score" if self.qual.nil?
 141     unless (Seq::SCORE_MIN .. Seq::SCORE_MAX).include? min_qual
 142       raise TrimError, "minimum quality value: #{min_qual} out of range #{Seq::SCORE_MIN} .. #{Seq::SCORE_MAX}"
 143     end
 144   end
 145
 146   # Inline C functions for speed below.
 147   inline do |builder|
 148     # Method for locating the right trim position and return this.
 149     builder.c %{
 150       VALUE trim_right_pos_c(
 151         VALUE _qual,        // quality score string
 152         VALUE _len,         // length of quality score string
 153         VALUE _min_qual,    // minimum quality score
 154         VALUE _min_len,     // minimum quality length
 155         VALUE _score_base   // score base
 156       )
 157       {
 158         unsigned char *qual       = StringValuePtr(_qual);
 159         unsigned int   len        = FIX2UINT(_len);
 160         unsigned int   min_qual   = FIX2UINT(_min_qual);
 161         unsigned int   min_len    = FIX2UINT(_min_len);
 162         unsigned int   score_base = FIX2UINT(_score_base);
 163
 164         unsigned int i = 0;
 165         unsigned int c = 0;
 166
 167         while (i < len)
 168         {
 169           c = 0;
 170
 171           while ((c < min_len) && ((c + i) < len) && (qual[len - (c + i) - 1] - score_base >= min_qual))
 172             c++;
 173
 174           if (c == min_len)
 175             return UINT2NUM(len - i);
 176           else
 177             i += c;
 178
 179           i++;
 180         }
 181
 182         return UINT2NUM(0);
 183       }
 184     }
 185
 186     # Method for locating the left trim position and return this.
 187     builder.c %{
 188       VALUE trim_left_pos_c(
 189         VALUE _qual,        // quality score string
 190         VALUE _len,         // length of quality score string
 191         VALUE _min_qual,    // minimum quality score
 192         VALUE _min_len,     // minimum quality length
 193         VALUE _score_base   // score base
 194       )
 195       {
 196         unsigned char *qual       = StringValuePtr(_qual);
 197         unsigned int   len        = FIX2UINT(_len);
 198         unsigned int   min_qual   = FIX2UINT(_min_qual);
 199         unsigned int   min_len    = FIX2UINT(_min_len);
 200         unsigned int   score_base = FIX2UINT(_score_base);
 201
 202         unsigned int i = 0;
 203         unsigned int c = 0;
 204
 205         while (i < len)
 206         {
 207           c = 0;
 208
 209           while ((c < min_len) && ((c + i) < len) && (qual[c + i] - score_base >= min_qual))
 210             c++;
 211
 212           if (c == min_len)
 213             return UINT2NUM(i);
 214           else
 215             i += c;
 216
 217           i++;
 218         }
 219
 220         return UINT2NUM(i);
 221       }
 222     }
 223   end
 224 end
 225
 226 __END__