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     @overlap = 0
  43     @offset1 = 0
  44     @offset2 = 0
  45     @options = options
  46     @options[:mismatches_max] ||= 0
  47     @options[:overlap_min]    ||= 1
  48   end
  49
  50   # Method to locate overlapping matches between two sequences.
  51   def match
  52     if @options[:overlap_max]
  53       @overlap = [@options[:overlap_max], @entry1.length, @entry2.length].min
  54     else
  55       @overlap = [@entry1.length, @entry2.length].min
  56     end
  57
  58     diff = @entry1.length - @entry2.length
  59     diff = 0 if diff < 0
  60
  61     @offset1 = @entry1.length - @overlap - diff
  62
  63     while @overlap >= @options[:overlap_min]
  64       mismatches_max = (@overlap * @options[:mismatches_max] * 0.01).round
  65
  66       # puts "diff: #{diff}   offset1: #{@offset1}  offset2: #{@offset2}   overlap: #{@overlap}"
  67
  68       if mismatches = match_C(@entry1.seq, @entry2.seq, @offset1, @offset2, @overlap, mismatches_max) and mismatches >= 0
  69         entry_merged          = entry_left + entry_overlap + entry_right
  70         entry_merged.seq_name = @entry1.seq_name + ":overlap=#{@overlap}:hamming=#{mismatches}"
  71
  72         return entry_merged
  73       end
  74
  75       if diff > 0
  76         diff -= 1
  77       else
  78         @overlap -= 1
  79       end
  80
  81       @offset1 += 1
  82     end
  83   end
  84
  85   # Method to extract and downcase the left part of an assembled pair.
  86   def entry_left
  87     entry = @entry1[0 ... @offset1]
  88     entry.seq.downcase!
  89     entry
  90   end
  91
  92   # Method to extract and downcase the right part of an assembled pair.
  93   def entry_right
  94     if @entry1.length > @offset1 + @overlap
  95       entry = @entry1[@offset1 + @overlap .. -1]
  96     else
  97       entry = @entry2[@offset2 + @overlap .. -1]
  98     end
  99
 100     entry.seq.downcase!
 101     entry
 102   end
 103
 104   # Method to extract and upcase the overlapping part of an assembled pair.
 105   def entry_overlap
 106     if @entry1.qual and @entry2.qual
 107       entry_overlap1 = @entry1[@offset1 ... @offset1 + @overlap]
 108       entry_overlap2 = @entry2[@offset2 ... @offset2 + @overlap]
 109
 110       entry = merge_overlap(entry_overlap1, entry_overlap2)
 111     else
 112       entry = @entry1[@offset1 ... @offset1 + @overlap]
 113     end
 114
 115     entry.seq.upcase!
 116     entry
 117   end
 118
 119   # Method to merge sequence and quality scores in an overlap.
 120   # The residue with the highest score at mismatch positions is selected.
 121   # The quality scores of the overlap are the mean of the two sequences.
 122   def merge_overlap(entry_overlap1, entry_overlap2)
 123     na_seq = NArray.byte(entry_overlap1.length, 2)
 124     na_seq[true, 0] = NArray.to_na(entry_overlap1.seq.downcase, "byte")
 125     na_seq[true, 1] = NArray.to_na(entry_overlap2.seq.downcase, "byte")
 126
 127     na_qual = NArray.byte(entry_overlap1.length, 2)
 128     na_qual[true, 0] = NArray.to_na(entry_overlap1.qual, "byte")
 129     na_qual[true, 1] = NArray.to_na(entry_overlap2.qual, "byte")
 130
 131     mask_xor = na_seq[true, 0] ^ na_seq[true, 1] > 0
 132     mask_seq = ((na_qual * mask_xor).eq( (na_qual * mask_xor).max(1)))
 133
 134     merged      = Seq.new()
 135     merged.seq  = (na_seq * mask_seq).max(1).to_s
 136     merged.qual = na_qual.mean(1).round.to_type("byte").to_s
 137
 138     merged
 139   end
 140
 141   inline do |builder|
 142     add_ambiguity_macro(builder)
 143
 144     # C method for determining if two strings of equal length match
 145     # given a maximum allowed mismatches and allowing for IUPAC
 146     # ambiguity codes. Returns number of mismatches is true if match, else false.
 147     builder.c %{
 148       VALUE match_C(
 149         VALUE _string1,       // String 1
 150         VALUE _string2,       // String 2
 151         VALUE _offset1,       // Offset 1
 152         VALUE _offset2,       // Offset 2
 153         VALUE _length,        // String length
 154         VALUE _max_mismatch   // Maximum mismatches
 155       )
 156       {
 157         char         *string1      = StringValuePtr(_string1);
 158         char         *string2      = StringValuePtr(_string2);
 159         unsigned int  offset1      = FIX2UINT(_offset1);
 160         unsigned int  offset2      = FIX2UINT(_offset2);
 161         unsigned int  length       = FIX2UINT(_length);
 162         unsigned int  max_mismatch = FIX2UINT(_max_mismatch);
 163
 164         unsigned int max_match = length - max_mismatch;
 165         unsigned int match     = 0;
 166         unsigned int mismatch  = 0;
 167         unsigned int i         = 0;
 168
 169         for (i = 0; i < length; i++)
 170         {
 171           if (MATCH(string1[i + offset1], string2[i + offset2]))
 172           {
 173             match++;
 174
 175             if (match >= max_match) {
 176               return UINT2NUM(mismatch);
 177             }
 178           }
 179           else
 180           {
 181             mismatch++;
 182
 183             if (mismatch > max_mismatch) {
 184               return INT2NUM(-1);
 185             }
 186           }
 187         }
 188
 189         return INT2NUM(-1);
 190       }
 191     }
 192   end
 193 end
 194
 195
 196 __END__
 197