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