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 # Module containing methods for end trimming sequences with suboptimal quality
  27 # scores.
  28 module Trim
  29   # Method to progressively trim a Seq object sequence from the right end until
  30   # a run of min_len residues with quality scores above min_qual is encountered.
  31   def quality_trim_right(min_qual, min_len = 1)
  32     check_trim_args(min_qual)
  33
  34     pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE)
  35
  36     self.subseq(0, pos)
  37   end
  38
  39   # Method to progressively trim a Seq object sequence from the right end until
  40   # a run of min_len residues with quality scores above min_qual is encountered.
  41   def quality_trim_right!(min_qual, min_len = 1)
  42     check_trim_args(min_qual)
  43
  44     pos = trim_right_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE)
  45
  46     self.subseq!(0, pos)
  47   end
  48
  49   # Method to progressively trim a Seq object sequence from the left end until
  50   # a run of min_len residues with quality scores above min_qual is encountered.
  51   def quality_trim_left(min_qual, min_len = 1)
  52     check_trim_args(min_qual)
  53
  54     pos = trim_left_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE)
  55
  56     self.subseq(pos)
  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, SCORE_BASE)
  65
  66     self.subseq!(pos)
  67   end
  68
  69   # Method to progressively trim a Seq object sequence from both ends until a
  70   # run of min_len residues with quality scores above min_qual is encountered.
  71   def quality_trim(min_qual, min_len = 1)
  72     check_trim_args(min_qual)
  73
  74     pos_right = trim_right_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE)
  75     pos_left  = trim_left_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE)
  76
  77     pos_left = pos_right if pos_left > pos_right
  78
  79     self.subseq(pos_left, pos_right - pos_left)
  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, SCORE_BASE)
  88     pos_left  = trim_left_pos_c(self.qual, self.length, min_qual, min_len, SCORE_BASE)
  89
  90     pos_left = pos_right if pos_left > pos_right
  91
  92     self.subseq!(pos_left, pos_right - pos_left)
  93   end
  94
  95   private
  96
  97   # Method to check the arguments for trimming and raise on bad sequence, qualities,
  98   # and min_qual.
  99   def check_trim_args(min_qual)
 100     raise SeqError, "no sequence"      if self.seq.nil?
 101     raise SeqError, "no quality score" if self.qual.nil?
 102     unless (SCORE_MIN .. SCORE_MAX).include? min_qual
 103       raise SeqError, "minimum quality value: #{min_qual} out of range #{SCORE_MIN} .. #{SCORE_MAX}"
 104     end
 105   end
 106
 107   # Inline C functions for speed below.
 108   inline do |builder|
 109     # Method for locating the right trim position and return this.
 110     builder.c %{
 111       VALUE trim_right_pos_c(
 112         VALUE _qual,        // quality score string
 113         VALUE _len,         // length of quality score string
 114         VALUE _min_qual,    // minimum quality score
 115         VALUE _min_len,     // minimum quality length
 116         VALUE _score_base   // score base
 117       )
 118       {
 119         unsigned char *qual       = StringValuePtr(_qual);
 120         unsigned int   len        = FIX2UINT(_len);
 121         unsigned int   min_qual   = FIX2UINT(_min_qual);
 122         unsigned int   min_len    = FIX2UINT(_min_len);
 123         unsigned int   score_base = FIX2UINT(_score_base);
 124
 125         unsigned int i = 0;
 126         unsigned int c = 0;
 127
 128         while (i < len)
 129         {
 130           c = 0;
 131
 132           while ((c < min_len) && ((c + i) < len) && (qual[len - (c + i) - 1] - score_base > min_qual))
 133             c++;
 134
 135           if (c == min_len)
 136             return UINT2NUM(len - i);
 137           else
 138             i += c;
 139
 140           i++;
 141         }
 142
 143         return UINT2NUM(0);
 144       }
 145     }
 146
 147     # Method for locating the left trim position and return this.
 148     builder.c %{
 149       VALUE trim_left_pos_c(
 150         VALUE _qual,        // quality score string
 151         VALUE _len,         // length of quality score string
 152         VALUE _min_qual,    // minimum quality score
 153         VALUE _min_len,     // minimum quality length
 154         VALUE _score_base   // score base
 155       )
 156       {
 157         unsigned char *qual       = StringValuePtr(_qual);
 158         unsigned int   len        = FIX2UINT(_len);
 159         unsigned int   min_qual   = FIX2UINT(_min_qual);
 160         unsigned int   min_len    = FIX2UINT(_min_len);
 161         unsigned int   score_base = FIX2UINT(_score_base);
 162
 163         unsigned int i = 0;
 164         unsigned int c = 0;
 165
 166         while (i < len)
 167         {
 168           c = 0;
 169
 170           while ((c < min_len) && ((c + i) < len) && (qual[c + i] - score_base > min_qual))
 171             c++;
 172
 173           if (c == min_len)
 174             return UINT2NUM(i);
 175           else
 176             i += c;
 177
 178           i++;
 179         }
 180
 181         return UINT2NUM(i);
 182       }
 183     }
 184   end
 185 end
 186
 187 __END__