1 # Copyright (C) 2007-2013 Martin A. Hansen.
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.
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.
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.
17 # http://www.gnu.org/copyleft/gpl.html
19 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
21 # This software is part of the Biopieces framework (www.biopieces.org).
23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
25 # Class containing methods to assemble two overlapping sequences into a single.
27 # Class method to assemble two Seq objects.
28 def self.pair(entry1, entry2, options = {})
29 assemble = self.new(entry1, entry2, options)
33 # Method to initialize an Assembly object.
34 def initialize(entry1, entry2, options)
38 @options[:mismatches_max] ||= 0
39 @options[:overlap_min] ||= 1
40 @options[:overlap_max] ||= entry1.length
41 @options[:overlap_max] = [@options[:overlap_max], entry1.length, entry2.length].min
44 # Method to locate overlapping matche between two sequences.
46 overlap = @options[:overlap_max]
48 na_seq1 = NArray.to_na(@entry1.seq.downcase, "byte")
49 na_seq2 = NArray.to_na(@entry2.seq.downcase, "byte")
51 while overlap >= @options[:overlap_min]
52 hamming_dist = (na_seq1[-1 * overlap .. -1] ^ na_seq2[0 ... overlap]).count_true
54 if hamming_dist <= (overlap * @options[:mismatches_max] * 0.01).round
55 entry_left = @entry1[0 ... @entry1.length - overlap]
56 entry_right = @entry2[overlap .. -1]
58 if @entry1.qual and @entry2.qual
59 entry_overlap1 = @entry1[-1 * overlap .. -1]
60 entry_overlap2 = @entry2[0 ... overlap]
62 entry_overlap = merge_overlap(entry_overlap1, entry_overlap2)
64 entry_overlap = @entry1[-1 * overlap .. -1]
67 entry_left.seq.downcase!
68 entry_overlap.seq.upcase!
69 entry_right.seq.downcase!
70 entry_merged = entry_left + entry_overlap + entry_right
71 entry_merged.seq_name = @entry1.seq_name + ":overlap=#{overlap}:hamming=#{hamming_dist}"
80 # Method to merge sequence and quality scores in an overlap.
81 # The residue with the highest score at mismatch positions is selected.
82 # The quality scores of the overlap are the mean of the two sequences.
83 def merge_overlap(entry_overlap1, entry_overlap2)
84 na_seq = NArray.byte(entry_overlap1.length, 2)
85 na_seq[true, 0] = NArray.to_na(entry_overlap1.seq.downcase, "byte")
86 na_seq[true, 1] = NArray.to_na(entry_overlap2.seq.downcase, "byte")
88 na_qual = NArray.byte(entry_overlap1.length, 2)
89 na_qual[true, 0] = NArray.to_na(entry_overlap1.qual, "byte")
90 na_qual[true, 1] = NArray.to_na(entry_overlap2.qual, "byte")
92 mask_xor = na_seq[true, 0] ^ na_seq[true, 1] > 0
93 mask_seq = ((na_qual * mask_xor).eq( (na_qual * mask_xor).max(1)))
96 merged.seq = (na_seq * mask_seq).max(1).to_s
97 merged.qual = na_qual.mean(1).round.to_type("byte").to_s