Class: ViralSeq::SeqHashPair

Inherits:
Object
  • Object
show all
Defined in:
lib/viral_seq/seq_hash_pair.rb

Overview

Class for paired-end sequences.

Examples:

initialize a new SeqHashPair object from a directory containing paired-end sequences

my_seqhashpair = ViralSeq::SeqHashPair.fa('my_seq_directory')

join the paired-end sequences with an overlap of 100 bp

my_seqhashpair.join1(100)

join the paired-end sequences with unknown overlap, each pair of sequences has its own overlap size

my_seqhashpair.join2(model: :indiv)

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(dna_hash = {}, title = "", file = []) ⇒ SeqHashPair

initialize SeqHashPair object with @dna_hash, @title and @file



16
17
18
19
20
 
# File 'lib/viral_seq/seq_hash_pair.rb', line 16

def initialize (dna_hash = {}, title = "", file = [])
  @dna_hash = dna_hash
  @title = title
  @file = file
end

Instance Attribute Details

#dna_hash ⇒ Hash

Returns Hash object for :name => [:r1_sequence_string, :r2_sequence_string].

Returns:

  • (Hash) —

    Hash object for :name => [:r1_sequence_string, :r2_sequence_string]



24
25
26
 
# File 'lib/viral_seq/seq_hash_pair.rb', line 24

def dna_hash
  @dna_hash
end

#file ⇒ String

Returns the r1 and r2 files that are used to initialize SeqHash object, if they exist.

Returns:

  • (String) —

    the r1 and r2 files that are used to initialize SeqHash object, if they exist



33
34
35
 
# File 'lib/viral_seq/seq_hash_pair.rb', line 33

def file
  @file
end

#title ⇒ String

default as the directory basename if SeqHash object is initialized using ::fa

Returns:

  • (String) —

    the title of the SeqHash object.



29
30
31
 
# File 'lib/viral_seq/seq_hash_pair.rb', line 29

def title
  @title
end

Class Method Details

.new_from_fasta(indir) ⇒ ViralSeq::SeqHashPair Also known as: fa

initialize a new ViralSeq::SeqHashPair object from a directory containing paired sequence files in the FASTA format

Examples:

initialize a new SeqHashPair object from a directory containing paired-end sequences

my_seqhashpair = ViralSeq::SeqHashPair.fa('spec/sample_paired_seq')

Parameters:

  • indir (String) —

    directory containing paired sequence files in the FASTA format,

    Paired sequence files need to have “r1” and “r2” in their file names

    Example for the file structure

    ├───lib1
    │     lib1_r1.txt
    │     lib1_r2.txt

    The sequence taxa should only differ by last 3 characters to distinguish r1 and r2 sequence.

Returns:



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
 
# File 'lib/viral_seq/seq_hash_pair.rb', line 49

def self.new_from_fasta(indir)
  files = Dir[indir + "/*"]
  r1_file = ""
  r2_file = ""
  files.each do |f|
    if File.basename(f) =~ /r1/i
      r1_file = f
    elsif File.basename(f) =~ /r2/i
      r2_file = f
    end
  end

  seq1 = ViralSeq::SeqHash.fa(r1_file).dna_hash
  seq2 = ViralSeq::SeqHash.fa(r2_file).dna_hash

  new_seq1 = seq1.each_with_object({}) {|(k, v), h| h[k[0..-4]] = v}
  new_seq2 = seq2.each_with_object({}) {|(k, v), h| h[k[0..-4]] = v}

  seq_pair_hash = {}

  new_seq1.each do |seq_name,seq|
    seq_pair_hash[seq_name] = [seq, new_seq2[seq_name]]
  end
  seq_hash = ViralSeq::SeqHashPair.new
  seq_hash.dna_hash = seq_pair_hash
  seq_hash.title = File.basename(indir,".*")
  seq_hash.file = [r1_file, r2_file]
  return seq_hash
end

Instance Method Details

#join1(overlap = 0, diff = 0.0) ⇒ ViralSeq::SeqHash

Pair-end join function for KNOWN overlap size. overlap can also be an explicit [Hash] object for :overlap_size, :r1_overlap, :r2_overlap, :before_overlap, :after_overlap

Examples:

join paired-end sequences with different :diff cut-offs, overlap provided.

paired_seqs = {">pair1"=>["GGGGGGGGGGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA",
                          "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATTTTTTTTTT"],
               ">pair2"=>["GGGGGGGGGGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA",
                          "AAAAAAAAAAGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATTTTTTTTTT"],
               ">pair3"=>["GGGGGGGGGGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA",
                          "AAAAAAAAAAGGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATTTTTTTTTT"]}
my_seqhashpair = ViralSeq::SeqHashPair.new(paired_seqs)
my_seqhashpair.join1(100).dna_hash.keys
=> [">pair1"]
my_seqhashpair.join1(100,0.01).dna_hash.keys
=> [">pair1", ">pair2"]
my_seqhashpair.join1(100,0.02).dna_hash.keys
=> [">pair1", ">pair2", ">pair3"]

Parameters:

  • overlap (Integer) (defaults to: 0) —

    simple overlap value indicating how many bases are overlapped. ‘0` means no overlap, R1 and R2 will be simply put together.

  • diff (Integer, Float) (defaults to: 0.0) —

    the maximum mismatch rate allowed for the overlapping region. default at 0.0, i.e. no mis-match allowed.

Returns:

Raises:

  • (ArgumentError) 


109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
 
# File 'lib/viral_seq/seq_hash_pair.rb', line 109

def join1(overlap = 0, diff = 0.0)
  raise ArgumentError.new(":diff has to be float or integer, input #{diff} invalid.") unless (diff.is_a? Integer or diff.is_a? Float)

  if overlap.is_a? Integer and overlap.zero?
    overlap = {
      overlap_size: 0,
      r1_overlapped: 0...0,
      r2_overlapped: 0...0,
      before_overlap: {
        region: :r1,
        range: 0..-1,
      } ,
      after_overlap: {
        region: :r2,
        range: 0..-1
      }
    }
  elsif overlap.is_a? Integer
    overlap = {
      overlap_size: overlap,
      r1_overlapped: -overlap..-1,
      r2_overlapped: 0..(overlap - 1),
      before_overlap: {
        region: :r1,
        range: 0..(-overlap - 1),
      } ,
      after_overlap: {
        region: :r2,
        range: overlap..-1
      }
    }
  end

  seq_pair_hash = self.dna_hash
  joined_seq = {}
  seq_pair_hash.each do |seq_name,seq_pair|
    r1_seq = seq_pair[0]
    r2_seq = seq_pair[1]

    r1_overlap = r1_seq[overlap[:r1_overlapped]]
    r2_overlap = r2_seq[overlap[:r2_overlapped]]

    overlap_size = overlap[:overlap_size]

    if (diff.zero? and r1_overlap == r2_overlap) or (!diff.zero? and r1_overlap.compare_with(r2_overlap) <= (overlap_size.abs * diff))
      if overlap[:before_overlap][:region] == :r1
        before_overlap_seq = r1_seq[overlap[:before_overlap][:range]]
      elsif overlap[:before_overlap][:region] == :r2
        before_overlap_seq = r2_seq[overlap[:before_overlap][:range]]
      end

      if overlap[:after_overlap][:region] == :r1
        after_overlap_seq = r1_seq[overlap[:after_overlap][:range]]
      elsif overlap[:after_overlap][:region] == :r2
        after_overlap_seq = r2_seq[overlap[:after_overlap][:range]]
      end
      joined_sequence = before_overlap_seq + r1_overlap + after_overlap_seq
    end

    joined_seq[seq_name] = joined_sequence if joined_sequence
  end

  joined_seq_hash = ViralSeq::SeqHash.new
  joined_seq_hash.dna_hash = joined_seq
  joined_seq_hash.title = self.title + "_joined"
  joined_seq_hash.file = File.dirname(self.file[0]) if self.file.size > 0
  return joined_seq_hash
end

#join2(model: :con, diff: 0.0) ⇒ ViralSeq::SeqHash

Pair-end join function for UNKNOWN overlap.

Examples:

join paired-end sequences, overlap NOT provided

paired_seq2 = {">pair4" => ["AAAGGGGGGG", "GGGGGGGTT"],
               ">pair5" => ["AAAAAAGGGG", "GGGGTTTTT"],
               ">pair6" => ["AAACAAGGGG", "GGGGTTTTT"] }
my_seqhashpair = ViralSeq::SeqHashPair.new(paired_seq2)
my_seqhashpair.join2.dna_hash
=> {">pair4"=>"AAAGGGGGGGGGGTT", ">pair5"=>"AAAAAAGGGGTTTTT", ">pair6"=>"AAACAAGGGGTTTTT"}
my_seqhashpair.join2(model: :indiv).dna_hash
=> {">pair4"=>"AAAGGGGGGGTT", ">pair5"=>"AAAAAAGGGGTTTTT", ">pair6"=>"AAACAAGGGGTTTTT"}

Parameters:

  • model (Symbol) (defaults to: :con) —

    models used to determine the overlap, ‘:con`, `:indiv`

    model ‘:con`: overlap is determined based on consensus, all sequence pairs are supposed to have the same overlap size

    note: minimal overlap as 4 bases.

    model ‘:indiv`: overlap is determined for each sequence pair, sequence pairs can have different size of overlap

  • diff (Integer, Float) (defaults to: 0.0) —

    the maximum mismatch rate allowed for the overlapping region. default at 0.0, i.e. no mis-match allowed.

Returns:



198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
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
247
248
249
 
# File 'lib/viral_seq/seq_hash_pair.rb', line 198

def join2(model: :con, diff: 0.0)
  seq_pair_hash = self.dna_hash
  begin
    raise ArgumentError.new(":diff has to be float or integer, input #{diff} invalid.") unless (diff.is_a? Integer or diff.is_a? Float)
    if model == :con
      overlap = determine_overlap_pid_pair(seq_pair_hash, diff)
      return self.join1(overlap, diff)
    elsif model == :indiv
      joined_seq = {}
      seq_pair_hash.each do |seq_name, seq_pair|
        r1_seq = seq_pair[0]
        r2_seq = seq_pair[1]
        overlap_list = []

        overlap_matrix(r1_seq, r2_seq).each do |overlap1, diff_nt|
          cut_off_base = overlap1[:overlap_size] * diff
          overlap_list << overlap1 if diff_nt <= cut_off_base
        end

        if overlap_list.empty?
          joined_seq[seq_name]  = seq_pair[0] + seq_pair[1]
        else
          overlap_to_use = overlap_list.sort_by{|k| k[:overlap_size].abs}.reverse[0]

          if overlap_to_use[:before_overlap][:region] == :r1
            before_overlap_seq = r1_seq[overlap_to_use[:before_overlap][:range]]
          elsif overlap_to_use[:before_overlap][:region] == :r2
            before_overlap_seq = r2_seq[overlap_to_use[:before_overlap][:range]]
          end

          if overlap_to_use[:after_overlap][:region] == :r1
            after_overlap_seq = r1_seq[overlap_to_use[:after_overlap][:range]]
          elsif overlap_to_use[:after_overlap][:region] == :r2
            after_overlap_seq = r2_seq[overlap_to_use[:after_overlap][:range]]
          end
          joined_seq[seq_name] = before_overlap_seq + r1_seq[overlap_to_use[:r1_overlapped]] + after_overlap_seq
        end
      end

      joined_seq_hash = ViralSeq::SeqHash.new
      joined_seq_hash.dna_hash = joined_seq
      joined_seq_hash.title = self.title + "_joined"
      joined_seq_hash.file = File.dirname(self.file[0]) if self.file.size > 0
      return joined_seq_hash
    else
      raise ArgumentError.new("Error::Wrong Overlap Model Argument. Given \`#{model}\`, expected `:con` or `:indiv`.")
    end
  rescue ArgumentError => e
    puts e
    return nil
  end
end

#size ⇒ Integer

the size of nt sequence hash of the SeqHashPair object

Returns:

  • (Integer) —

    size of nt sequence hash of the SeqHash object



85
86
87
 
# File 'lib/viral_seq/seq_hash_pair.rb', line 85

def size
  self.dna_hash.size
end