code_ruby/lib/maasha/seq/trim.rb

   1
   2 # This program is free software; you can redistribute it and/or
   3 # modify it under the terms of the GNU General Public License
   4 # as published by the Free Software Foundation; either version 2
   5 # of the License, or (at your option) any later version.
   6
   7 # This program is distributed in the hope that it will be useful,
   8 # but WITHOUT ANY WARRANTY; without even the implied warranty of
   9 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  10 # GNU General Public License for more details.
  11
  12 # You should have received a copy of the GNU General Public License
  13 # along with this program; if not, write to the Free Software
  14 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  15
  16 # http://www.gnu.org/copyleft/gpl.html
  17
  18 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  19
  20 # This software is part of the Biopieces framework (www.biopieces.org).
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24 require 'inline'
  25
  26 # Error class for all exceptions to do with Trim.
  27 class TrimError < StandardError; end
  28
  29 # Module containing methods for end trimming sequences with suboptimal quality
  30 # scores.
  31 module Trim
  32   # Method to progressively trim a Seq object sequence from the right end until
  33   # a run of min_len residues with quality scores above min_qual is encountered.
  34   def quality_trim_right(min_qual, min_len = 1)
  35     check_trim_args(min_qual)
  36
  37     pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  38
  39     self.subseq(0, pos)
  40   end
  41
  42   # Method to progressively trim a Seq object sequence from the right end until
  43   # a run of min_len residues with quality scores above min_qual is encountered.
  44   def quality_trim_right!(min_qual, min_len = 1)
  45     check_trim_args(min_qual)
  46
  47     pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  48
  49     self.subseq!(0, pos)
  50   end
  51
  52   # Method to progressively trim a Seq object sequence from the left end until
  53   # a run of min_len residues with quality scores above min_qual is encountered.
  54   def quality_trim_left(min_qual, min_len = 1)
  55     check_trim_args(min_qual)
  56
  57     pos = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  58
  59     self.subseq(pos)
  60   end
  61
  62   # Method to progressively trim a Seq object sequence from the left end until
  63   # a run of min_len residues with quality scores above min_qual is encountered.
  64   def quality_trim_left!(min_qual, min_len = 1)
  65     check_trim_args(min_qual)
  66
  67     pos = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  68
  69     self.subseq!(pos)
  70   end
  71
  72   # Method to progressively trim a Seq object sequence from both ends until a
  73   # run of min_len residues with quality scores above min_qual is encountered.
  74   def quality_trim(min_qual, min_len = 1)
  75     check_trim_args(min_qual)
  76
  77     pos_right = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  78     pos_left  = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  79
  80     pos_left = pos_right if pos_left > pos_right
  81
  82     self.subseq(pos_left, pos_right - pos_left)
  83   end
  84
  85   # Method to progressively trim a Seq object sequence from both ends until a
  86   # run of min_len residues with quality scores above min_qual is encountered.
  87   def quality_trim!(min_qual, min_len = 1)
  88     check_trim_args(min_qual)
  89
  90     pos_right = trim_right_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  91     pos_left  = trim_left_pos_c(self.qual, self.length, min_qual, min_len, Seq::SCORE_BASE)
  92
  93     pos_left = pos_right if pos_left > pos_right
  94
  95     self.subseq!(pos_left, pos_right - pos_left)
  96   end
  97
  98   private
  99
 100   # Method to check the arguments for trimming and raise on bad sequence, qualities,
 101   # and min_qual.
 102   def check_trim_args(min_qual)
 103     raise TrimError, "no sequence"      if self.seq.nil?
 104     raise TrimError, "no quality score" if self.qual.nil?
 105     unless (Seq::SCORE_MIN .. Seq::SCORE_MAX).include? min_qual
 106       raise TrimError, "minimum quality value: #{min_qual} out of range #{Seq::SCORE_MIN} .. #{Seq::SCORE_MAX}"
 107     end
 108   end
 109
 110   # Inline C functions for speed below.
 111   inline do |builder|
 112     # Method for locating the right trim position and return this.
 113     builder.c %{
 114       VALUE trim_right_pos_c(
 115         VALUE _qual,        // quality score string
 116         VALUE _len,         // length of quality score string
 117         VALUE _min_qual,    // minimum quality score
 118         VALUE _min_len,     // minimum quality length
 119         VALUE _score_base   // score base
 120       )
 121       {
 122         unsigned char *qual       = StringValuePtr(_qual);
 123         unsigned int   len        = FIX2UINT(_len);
 124         unsigned int   min_qual   = FIX2UINT(_min_qual);
 125         unsigned int   min_len    = FIX2UINT(_min_len);
 126         unsigned int   score_base = FIX2UINT(_score_base);
 127
 128         unsigned int i = 0;
 129         unsigned int c = 0;
 130
 131         while (i < len)
 132         {
 133           c = 0;
 134
 135           while ((c < min_len) && ((c + i) < len) && (qual[len - (c + i) - 1] - score_base >= min_qual))
 136             c++;
 137
 138           if (c == min_len)
 139             return UINT2NUM(len - i);
 140           else
 141             i += c;
 142
 143           i++;
 144         }
 145
 146         return UINT2NUM(0);
 147       }
 148     }
 149
 150     # Method for locating the left trim position and return this.
 151     builder.c %{
 152       VALUE trim_left_pos_c(
 153         VALUE _qual,        // quality score string
 154         VALUE _len,         // length of quality score string
 155         VALUE _min_qual,    // minimum quality score
 156         VALUE _min_len,     // minimum quality length
 157         VALUE _score_base   // score base
 158       )
 159       {
 160         unsigned char *qual       = StringValuePtr(_qual);
 161         unsigned int   len        = FIX2UINT(_len);
 162         unsigned int   min_qual   = FIX2UINT(_min_qual);
 163         unsigned int   min_len    = FIX2UINT(_min_len);
 164         unsigned int   score_base = FIX2UINT(_score_base);
 165
 166         unsigned int i = 0;
 167         unsigned int c = 0;
 168
 169         while (i < len)
 170         {
 171           c = 0;
 172
 173           while ((c < min_len) && ((c + i) < len) && (qual[c + i] - score_base >= min_qual))
 174             c++;
 175
 176           if (c == min_len)
 177             return UINT2NUM(i);
 178           else
 179             i += c;
 180
 181           i++;
 182         }
 183
 184         return UINT2NUM(i);
 185       }
 186     }
 187   end
 188 end
 189
 190 __END__