code_ruby/lib/maasha/seq/assemble.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/ambiguity'
  27
  28 # Class containing methods to assemble two overlapping sequences into a single.
  29 class Assemble
  30   extend Ambiguity
  31
  32   # Class method to assemble two Seq objects.
  33   def self.pair(entry1, entry2, options = {})
  34     assemble = self.new(entry1, entry2, options)
  35     assemble.match
  36   end
  37
  38   # Method to initialize an Assembly object.
  39   def initialize(entry1, entry2, options)
  40     @entry1  = entry1
  41     @entry2  = entry2
  42     @options = options
  43     @options[:mismatches_max] ||= 0
  44     @options[:overlap_min]    ||= 1
  45     @options[:overlap_max]    ||= entry1.length
  46     @options[:overlap_max]      = [@options[:overlap_max], entry1.length, entry2.length].min
  47   end
  48
  49   # Method to locate overlapping matches between two sequences.
  50   def match
  51     overlap = @options[:overlap_max]
  52
  53     while overlap >= @options[:overlap_min]
  54       mismatches_max = (overlap * @options[:mismatches_max] * 0.01).round
  55
  56       if match_C(@entry1.seq, @entry2.seq, @entry1.length - overlap, 0, overlap, mismatches_max)
  57         entry_left  = @entry1[0 ... @entry1.length - overlap]
  58         entry_right = @entry2[overlap .. -1]
  59
  60         if @entry1.qual and @entry2.qual
  61           entry_overlap1 = @entry1[-1 * overlap .. -1]
  62           entry_overlap2 = @entry2[0 ... overlap]
  63
  64           entry_overlap = merge_overlap(entry_overlap1, entry_overlap2)
  65         else
  66           entry_overlap = @entry1[-1 * overlap .. -1]
  67         end
  68
  69         entry_left.seq.downcase!
  70         entry_overlap.seq.upcase!
  71         entry_right.seq.downcase!
  72         entry_merged          = entry_left + entry_overlap + entry_right
  73         entry_merged.seq_name = @entry1.seq_name + ":overlap=#{overlap}"
  74
  75         return entry_merged
  76       end
  77
  78       overlap -= 1
  79     end
  80   end
  81
  82   # Method to merge sequence and quality scores in an overlap.
  83   # The residue with the highest score at mismatch positions is selected.
  84   # The quality scores of the overlap are the mean of the two sequences.
  85   def merge_overlap(entry_overlap1, entry_overlap2)
  86     na_seq = NArray.byte(entry_overlap1.length, 2)
  87     na_seq[true, 0] = NArray.to_na(entry_overlap1.seq.downcase, "byte")
  88     na_seq[true, 1] = NArray.to_na(entry_overlap2.seq.downcase, "byte")
  89
  90     na_qual = NArray.byte(entry_overlap1.length, 2)
  91     na_qual[true, 0] = NArray.to_na(entry_overlap1.qual, "byte")
  92     na_qual[true, 1] = NArray.to_na(entry_overlap2.qual, "byte")
  93
  94     mask_xor = na_seq[true, 0] ^ na_seq[true, 1] > 0
  95     mask_seq = ((na_qual * mask_xor).eq( (na_qual * mask_xor).max(1)))
  96
  97     merged      = Seq.new()
  98     merged.seq  = (na_seq * mask_seq).max(1).to_s
  99     merged.qual = na_qual.mean(1).round.to_type("byte").to_s
 100
 101     merged
 102   end
 103
 104   inline do |builder|
 105     add_ambiguity_macro(builder)
 106
 107     # C method for determining if two strings of equal length match
 108     # given a maximum allowed mismatches and allowing for IUPAC
 109     # ambiguity codes. Returns true if match, else false.
 110     builder.c %{
 111       VALUE match_C(
 112         VALUE _string1,       // String 1
 113         VALUE _string2,       // String 2
 114         VALUE _offset1,       // Offset 1
 115         VALUE _offset2,       // Offset 2
 116         VALUE _length,        // String length
 117         VALUE _max_mismatch   // Maximum mismatches
 118       )
 119       {
 120         char         *string1      = StringValuePtr(_string1);
 121         char         *string2      = StringValuePtr(_string2);
 122         unsigned int  offset1      = FIX2UINT(_offset1);
 123         unsigned int  offset2      = FIX2UINT(_offset2);
 124         unsigned int  length       = FIX2UINT(_length);
 125         unsigned int  max_mismatch = FIX2UINT(_max_mismatch);
 126
 127         unsigned int max_match = length - max_mismatch;
 128         unsigned int match     = 0;
 129         unsigned int mismatch  = 0;
 130         unsigned int i         = 0;
 131
 132         for (i = 0; i < length; i++)
 133         {
 134           if (MATCH(string1[i + offset1], string2[i + offset2]))
 135           {
 136             match++;
 137
 138             if (match >= max_match) {
 139               return Qtrue;
 140             }
 141           }
 142           else
 143           {
 144             mismatch++;
 145
 146             if (mismatch > max_mismatch) {
 147               return Qfalse;
 148             }
 149           }
 150         }
 151
 152         return Qfalse;
 153       }
 154     }
 155   end
 156 end
 157
 158
 159 __END__
 160