Module: SSW

Defined in:
lib/libssw.rb,
lib/ssw/dna.rb,
lib/ssw/ffi.rb,
lib/ssw/aaseq.rb,
lib/ssw/align.rb,
lib/ssw/profile.rb,
lib/ssw/version.rb,
lib/ssw/BLOSUM50.rb,
lib/ssw/BLOSUM62.rb

Defined Under Namespace

Modules: AASeq, DNA, FFI Classes: Align, Error, Profile

Constant Summary collapse

VERSION = 
'0.0.3'
BLOSUM50 = 
[
  #  A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V   B   Z   X   *
 5, -2, -1, -2, -1, -1, -1,  0, -2, -1, -2, -1, -1, -3, -1,  1,  0, -3, -2,  0, -2, -1, -1, -5, # A
-2,  7, -1, -2, -4,  1,  0, -3,  0, -4, -3,  3, -2, -3, -3, -1, -1, -3, -1, -3, -1,  0, -1, -5, # R
-1, -1,  7,  2, -2,  0,  0,  0,  1, -3, -4,  0, -2, -4, -2,  1,  0, -4, -2, -3,  5,  0, -1, -5, # N
-2, -2,  2,  8, -4,  0,  2, -1, -1, -4, -4, -1, -4, -5, -1,  0, -1, -5, -3, -4,  6,  1, -1, -5, # D
-1, -4, -2, -4, 13, -3, -3, -3, -3, -2, -2, -3, -2, -2, -4, -1, -1, -5, -3, -1, -3, -3, -1, -5, # C
-1,  1,  0,  0, -3,  7,  2, -2,  1, -3, -2,  2,  0, -4, -1,  0, -1, -1, -1, -3,  0,  4, -1, -5, # Q
-1,  0,  0,  2, -3,  2,  6, -3,  0, -4, -3,  1, -2, -3, -1, -1, -1, -3, -2, -3,  1,  5, -1, -5, # E
 0, -3,  0, -1, -3, -2, -3,  8, -2, -4, -4, -2, -3, -4, -2,  0, -2, -3, -3, -4, -1, -2, -1, -5, # G
-2,  0,  1, -1, -3,  1,  0, -2, 10, -4, -3,  0, -1, -1, -2, -1, -2, -3,  2, -4,  0,  0, -1, -5, # H
-1, -4, -3, -4, -2, -3, -4, -4, -4,  5,  2, -3,  2,  0, -3, -3, -1, -3, -1,  4, -4, -3, -1, -5, # I
-2, -3, -4, -4, -2, -2, -3, -4, -3,  2,  5, -3,  3,  1, -4, -3, -1, -2, -1,  1, -4, -3, -1, -5, # L
-1,  3,  0, -1, -3,  2,  1, -2,  0, -3, -3,  6, -2, -4, -1,  0, -1, -3, -2, -3,  0,  1, -1, -5, # K
-1, -2, -2, -4, -2,  0, -2, -3, -1,  2,  3, -2,  7,  0, -3, -2, -1, -1,  0,  1, -3, -1, -1, -5, # M
-3, -3, -4, -5, -2, -4, -3, -4, -1,  0,  1, -4,  0,  8, -4, -3, -2,  1,  4, -1, -4, -4, -1, -5, # F
-1, -3, -2, -1, -4, -1, -1, -2, -2, -3, -4, -1, -3, -4, 10, -1, -1, -4, -3, -3, -2, -1, -1, -5, # P
 1, -1,  1,  0, -1,  0, -1,  0, -1, -3, -3,  0, -2, -3, -1,  5,  2, -4, -2, -2,  0,  0, -1, -5, # S
 0, -1,  0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1,  2,  5, -3, -2,  0,  0, -1, -1, -5, # T
-3, -3, -4, -5, -5, -1, -3, -3, -3, -3, -2, -3, -1,  1, -4, -4, -3, 15,  2, -3, -5, -2, -1, -5, # W
-2, -1, -2, -3, -3, -1, -2, -3,  2, -1, -1, -2,  0,  4, -3, -2, -2,  2,  8, -1, -3, -2, -1, -5, # Y
 0, -3, -3, -4, -1, -3, -3, -4, -4,  4,  1, -3,  1, -1, -3, -2,  0, -3, -1,  5, -3, -3, -1, -5, # V
-2, -1,  5,  6, -3,  0,  1, -1,  0, -4, -4,  0, -3, -4, -2,  0,  0, -5, -3, -3,  6,  1, -1, -5, # B
-1,  0,  0,  1, -3,  4,  5, -2,  0, -3, -3,  1, -1, -4, -1,  0, -1, -2, -2, -3,  1,  5, -1, -5, # Z
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -5, # X
-5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5,  1  # *
]
BLOSUM62 = 
[
  #  A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V   B   Z   X   *
 4, -1, -2, -2,  0, -1, -1,  0, -2, -1, -1, -1, -1, -2, -1,  1,  0, -3, -2,  0, -2, -1,  0, -4, # A
-1,  5,  0, -2, -3,  1,  0, -2,  0, -3, -2,  2, -1, -3, -2, -1, -1, -3, -2, -3, -1,  0, -1, -4, # R
-2,  0,  6,  1, -3,  0,  0,  0,  1, -3, -3,  0, -2, -3, -2,  1,  0, -4, -2, -3,  3,  0, -1, -4, # N
-2, -2,  1,  6, -3,  0,  2, -1, -1, -3, -4, -1, -3, -3, -1,  0, -1, -4, -3, -3,  4,  1, -1, -4, # D
 0, -3, -3, -3,  9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1, -2, -2, -1, -3, -3, -2, -4, # C
-1,  1,  0,  0, -3,  5,  2, -2,  0, -3, -2,  1,  0, -3, -1,  0, -1, -2, -1, -2,  0,  3, -1, -4, # Q
-1,  0,  0,  2, -4,  2,  5, -2,  0, -3, -3,  1, -2, -3, -1,  0, -1, -3, -2, -2,  1,  4, -1, -4, # E
 0, -2,  0, -1, -3, -2, -2,  6, -2, -4, -4, -2, -3, -3, -2,  0, -2, -2, -3, -3, -1, -2, -1, -4, # G
-2,  0,  1, -1, -3,  0,  0, -2,  8, -3, -3, -1, -2, -1, -2, -1, -2, -2,  2, -3,  0,  0, -1, -4, # H
-1, -3, -3, -3, -1, -3, -3, -4, -3,  4,  2, -3,  1,  0, -3, -2, -1, -3, -1,  3, -3, -3, -1, -4, # I
-1, -2, -3, -4, -1, -2, -3, -4, -3,  2,  4, -2,  2,  0, -3, -2, -1, -2, -1,  1, -4, -3, -1, -4, # L
-1,  2,  0, -1, -3,  1,  1, -2, -1, -3, -2,  5, -1, -3, -1,  0, -1, -3, -2, -2,  0,  1, -1, -4, # K
-1, -1, -2, -3, -1,  0, -2, -3, -2,  1,  2, -1,  5,  0, -2, -1, -1, -1, -1,  1, -3, -1, -1, -4, # M
-2, -3, -3, -3, -2, -3, -3, -3, -1,  0,  0, -3,  0,  6, -4, -2, -2,  1,  3, -1, -3, -3, -1, -4, # F
-1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4,  7, -1, -1, -4, -3, -2, -2, -1, -2, -4, # P
 1, -1,  1,  0, -1,  0,  0,  0, -1, -2, -2,  0, -1, -2, -1,  4,  1, -3, -2, -2,  0,  0,  0, -4, # S
 0, -1,  0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1,  1,  5, -2, -2,  0, -1, -1,  0, -4, # T
-3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1,  1, -4, -3, -2, 11,  2, -3, -4, -3, -2, -4, # W
-2, -2, -2, -3, -2, -1, -2, -3,  2, -1, -1, -2, -1,  3, -3, -2, -2,  2,  7, -1, -3, -2, -1, -4, # Y
 0, -3, -3, -3, -1, -2, -2, -3, -3,  3,  1, -2,  1, -1, -2, -2,  0, -3, -1,  4, -3, -2, -1, -4, # V
-2, -1,  3,  4, -3,  0,  1, -1,  0, -3, -4,  0, -3, -3, -2,  0, -1, -4, -3, -3,  4,  1, -1, -4, # B
-1,  0,  0,  1, -3,  3,  4, -2,  0, -3, -3,  1, -1, -3, -1,  0, -1, -3, -2, -2,  1,  4, -1, -4, # Z
 0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2,  0,  0, -2, -1, -1, -1, -1, -1, -4, # X
-4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4,  1, # *
]

Class Attribute Summary collapse

Class Method Summary collapse

Class Attribute Details

.ffi_lib ⇒ Object

Returns the value of attribute ffi_lib.



13
14
15
 
# File 'lib/libssw.rb', line 13

def ffi_lib
  @ffi_lib
end

Class Method Details

.align(prof, ref, weight_gap0, weight_gapE, flag, filters, filterd, mask_len = nil) ⇒ Object

Do Striped Smith-Waterman alignment.

Parameters:

  • prof (Fiddle::Pointer, SSW::Profile, SSW::FFI::Profile) —

    pointer to the query profile structure

  • ref (Array) —

    target sequence; the target sequence needs to be numbers and corresponding to the mat parameter of function ssw_init

  • weight_gap0 (Integer) —

    the absolute value of gap open penalty

  • weight_gapE (Integer) —

    the absolute value of gap extension penalty

  • flag (Integer) —

    • bit 5: when setted as 1, function ssw_align will return the best alignment beginning position;

    • bit 6: when setted as 1, if (ref_end1 - ref_begin1 < filterd && read_end1 - read_begin1 < filterd), (whatever bit 5 is setted) the function will return the best alignment beginning position and cigar;

    • bit 7: when setted as 1, if the best alignment score >= filters, (whatever bit 5 is setted) the function will return the best alignment beginning position and cigar;

    • bit 8: when setted as 1, (whatever bit 5, 6 or 7 is setted) the function will always return the best alignment beginning position and cigar. When flag == 0, only the optimal and sub-optimal scores and the optimal alignment ending position will be returned.

  • filters (Integer) —

    scorefilter: when bit 7 of flag is setted as 1 and bit 8 is setted as 0, filters will be used (Please check the decription of the flag parameter for detailed usage.)

  • filterd (Integer) —

    distance filter: when bit 6 of flag is setted as 1 and bit 8 is setted as 0, filterd will be used (Please check the decription of the flag parameter for detailed usage.)

  • mask_len (Integer) (defaults to: nil) —

    The distance between the optimal and suboptimal alignment ending position >= maskLen. We suggest to use readLen/2, if you don’t have special concerns. Note: maskLen has to be >= 15, otherwise this function will NOT return the suboptimal alignment information. Detailed description of maskLen: After locating the optimal alignment ending position, the suboptimal alignment score can be heuristically found by checking the second largest score in the array that contains the maximal score of each column of the SW matrix. In order to avoid picking the scores that belong to the alignments sharing the partial best alignment, SSW C library masks the reference loci nearby (mask length = maskLen) the best alignment ending position and locates the second largest score from the unmasked elements.

Raises:

  • (ArgumentError) 


136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
 
# File 'lib/libssw.rb', line 136

def align(prof, ref, weight_gap0, weight_gapE, flag, filters, filterd, mask_len = nil)
  unless prof.is_a?(Fiddle::Pointer) || prof.is_a?(Profile) || prof.respond_to?(:to_ptr)
    raise ArgumentError, 'Expect class of filename to be Profile or Pointer'
  end
  raise ArgumentError, 'Expect class of ref to be Array' unless ref.is_a?(Array)

  ref_str = ref.pack('c*')
  ref_len = ref.size
  mask_len ||= [ref_len / 2, 15].max
  ptr = FFI.ssw_align(
    prof, ref_str, ref_len, weight_gap0, weight_gapE, flag, filters, filterd, mask_len
  )
  # Not sure yet if we should set the instance variable to the pointer as a
  # garbage collection workaround.
  # For example: instance_variable_set(:@ref_str, ref_str)
  #
  # ptr.free = FFI.instance_variable_get(:@func_map)['align_destroy']
  SSW::Align.new(ptr)
end

.align_destroy(align) ⇒ Object

Release the memory allocated by function ssw_align.

Parameters:



159
160
161
162
163
164
165
166
167
168
 
# File 'lib/libssw.rb', line 159

def align_destroy(align)
  if align.is_a?(Align)
    warn "You don't need to call this method for Ruby's Align class."
    nil
  elsif align.is_a?(Fiddle::Pointer) || align.respond_to?(:to_ptr)
    FFI.align_destroy(align)
  else
    raise ArgumentError, 'Expect class of filename to be Pointer'
  end
end

.build_path(q_seq, r_seq, align) ⇒ Object

TODO: fix variable names



216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
 
# File 'lib/libssw.rb', line 216

def build_path(q_seq, r_seq, align)
  sQ = ''
  sA = ''
  sR = ''
  q_off = align.read_begin1
  r_off = align.ref_begin1
  align.cigar.each do |x|
    n = x >> 4
    m = x & 15
    c = m > 8 ? 'M' : 'MIDNSHP=X'[m]
    case c
    when 'M'
      sQ += q_seq[q_off...(q_off + n)]
      sA += Array.new(n) { |j| q_seq[q_off + j] == r_seq[r_off + j] ? '|' : '*' }.join
      sR += r_seq[r_off...(r_off + n)]
      q_off += n
      r_off += n
    when 'I'
      sQ += q_seq[q_off...(q_off + n)]
      sA += ' ' * n
      sR += ' ' * n
      q_off += n
    when 'D'
      sQ += ' ' * n
      sA += ' ' * n
      sR += r_seq[r_off...(r_off + n)]
      r_off += n
    end
  end
  [align.cigar_string, sQ, sA, sR]
end

.create_scoring_matrix(elements, match_score, mismatch_score) ⇒ Object

Create scoring matrix of Smith-Waterman algrithum.

Parameters:

  • elements (Array)
  • match_score (Integer)
  • mismatch_score (Integer) 


203
204
205
206
207
208
209
210
211
212
213
 
# File 'lib/libssw.rb', line 203

def create_scoring_matrix(elements, match_score, mismatch_score)
  size = elements.size
  score = Array.new(size * size, 0)
  (size - 1).times do |i|
    (size - 1).times do |j|
      score[i * size + j] = \
        (elements[i] == elements[j] ? match_score : mismatch_score)
    end
  end
  score
end

.init(read, mat, n = nil, score_size: 2) ⇒ Object

Create the query profile using the query sequence.

Parameters:

  • read (Array) —

    query sequence; the query sequence needs to be numbers

  • mat (Array) —

    substitution matrix; mat needs to be corresponding to the read sequence

  • n (Integer) (defaults to: nil) —

    the square root of the number of elements in mat (mat has n*n elements) If you omit this argument, the square root of the size of mat will be set.

  • score_size (Integer) (defaults to: 2) —

    estimated Smith-Waterman score;

    • if your estimated best alignment score is surely < 255 please set 0;

    • if your estimated best alignment score >= 255, please set 1;

    • if you don’t know, please set 2

Raises:

  • (ArgumentError) 


46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
 
# File 'lib/libssw.rb', line 46

def init(read, mat, n = nil, score_size: 2)
  read = read.to_a
  mat = mat.to_a.flatten
  raise ArgumentError, 'Expect class of read to be Array' unless read.is_a?(Array)
  raise ArgumentError, 'Expect class of mat to be Array' unless mat.is_a?(Array)

  read_str = read.pack('c*')
  read_len = read.size
  n = Math.sqrt(mat.size) if n.nil?
  raise "Not a square matrix. size: #{mat.size}, n: #{n}" if mat.size != n * n

  mat_str = mat.flatten.pack('c*')
  ptr = FFI.ssw_init(
    read_str,
    read_len,
    mat_str,
    n,
    score_size
  )
  # Garbage collection workaround
  #
  # * The following code will cause a segmentation violation when manually
  #   releasing memory. The reason is unknown.
  # * func_map is only available in newer versions of fiddle.
  # ptr.free = FFI.instance_variable_get(:@func_map)['init_destroy']
  ptr.instance_variable_set(:@read_str,   read_str)
  ptr.instance_variable_set(:@read_len,   read_len)
  ptr.instance_variable_set(:@mat_str,    mat_str)
  ptr.instance_variable_set(:@n,          n)
  ptr.instance_variable_set(:@score_size, score_size)

  SSW::Profile.new(ptr)
end

.init_destroy(profile) ⇒ Object

Note:

Ruby has garbage collection, so there is not much reason to call this method.

Release the memory allocated by function ssw_init.

Parameters:



85
86
87
88
89
90
91
 
# File 'lib/libssw.rb', line 85

def init_destroy(profile)
  unless profile.is_a?(Fiddle::Pointer) || prof.is_a?(Profile) || prof.respond_to?(:to_ptr)
    raise ArgumentError, 'Expect class of filename to be Profile or Pointer'
  end

  FFI.init_destroy(profile)
end

.mark_mismatch(ref_begin1, read_begin1, read_end1, ref, read, read_len, cigar, cigar_len) ⇒ Integer

Note:

This method takes a Fiddle::Pointer as an argument. Please read the source code and understand it well before using this method. (Needs to be improved)

  1. Calculate the number of mismatches.

  2. Modify the cigar string:

differentiate matches (=), mismatches(X), and softclip(S).

Parameters:

  • ref_begin1 (Integer) —

    0-based best alignment beginning position on the reference sequence

  • read_begin1 (Integer) —

    0-based best alignment beginning position on the read sequence

  • read_end1 (Integer) —

    0-based best alignment ending position on the read sequence

  • ref (Array) —

    reference sequence

  • read (Array) —

    read sequence

  • read_len (Integer) —

    length of the read

  • cigar (Fiddle::Pointer) —

    best alignment cigar; stored the same as that in BAM format, high 28 bits: length, low 4 bits: M/I/D (0/1/2)

  • cigar_len (Integer) —

    length of the cigar string

Returns:

  • (Integer) —

    The number of mismatches. The cigar and cigarLen are modified.



192
193
194
195
196
197
 
# File 'lib/libssw.rb', line 192

def mark_mismatch(ref_begin1, read_begin1, read_end1, ref, read, read_len, cigar, cigar_len)
  warn 'implementation: fiexme: **cigar' # FIXME
  FFI.mark_mismatch(
    ref_begin1, read_begin1, read_end1, ref.pack('c*'), read.pack('c*'), read_len, cigar, cigar_len.pack('l*')
  )
end