Class: Genfrag::App::SearchCommand

Inherits:
Command
  • Object
show all
Defined in:
lib/genfrag/app/search_command.rb,
lib/genfrag/app/search_command/trim.rb,
lib/genfrag/app/search_command/match.rb,
lib/genfrag/app/search_command/process_file.rb

Defined Under Namespace

Classes: ProcessFile

Instance Attribute Summary

Attributes inherited from Command

#ops, #options

Instance Method Summary collapse

Methods inherited from Command

#cli_p, #clierr_p, #initialize, #standard_options

Constructor Details

This class inherits a constructor from Genfrag::App::Command

Instance Method Details

#adapter_setup_1(hsh) ⇒ Object 



299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
 
# File 'lib/genfrag/app/search_command.rb', line 299

def adapter_setup_1(hsh)
  l = lambda do |i|
    if @ops.send("adapter#{i}")
      @adapters["adapter#{i}_specificity".to_sym] = @ops.send("adapter#{i}")
      if @ops.send("adapter#{i}_sequence")
        @adapters["adapter#{i}_sequence".to_sym] = @ops.send("adapter#{i}_sequence").gsub(/\|N*$/i,'')
        @adapters["adapter#{i}_size".to_sym] = @adapters["adapter#{i}_sequence".to_sym].size + @adapters["adapter#{i}_specificity".to_sym].size
      else
        @adapters["adapter#{i}_size".to_sym] = @ops.send("adapter#{i}_size")
      end
    elsif hsh["adapter#{i}_specificity".to_sym]
      @adapters["adapter#{i}_specificity".to_sym] = hsh["adapter#{i}_specificity".to_sym]
      @adapters["adapter#{i}_sequence".to_sym]    = hsh["adapter#{i}_sequence".to_sym]
      @adapters["adapter#{i}_size".to_sym]        = hsh["adapter#{i}_sequence".to_sym].size + hsh["adapter#{i}_specificity".to_sym].size
    end
  end
# set adapter 5' and 3' respectively using above procs
  l.call(5)
  l.call(3)
end

#adapter_setup_2 ⇒ Object 



320
321
322
323
324
325
326
327
328
329
330
331
332
 
# File 'lib/genfrag/app/search_command.rb', line 320

def adapter_setup_2
  l = lambda do |i|
    @adapters["adapter#{i}_specificity".to_sym] = @ops.send("adapter#{i}")
    if @ops.send("adapter#{i}_sequence")
      @adapters["adapter#{i}_sequence".to_sym] = @ops.send("adapter#{i}_sequence").gsub(/\|N*$/i,'')
      @adapters["adapter#{i}_size".to_sym] = @adapters["adapter#{i}_sequence".to_sym].size + @adapters["adapter#{i}_specificity".to_sym].size
    else
      @adapters["adapter#{i}_size".to_sym] = @ops.send("adapter#{i}_size")
    end
  end
  l.call(5)
  l.call(3)
end

#calculate_left_and_right_trims(trim) ⇒ Object

Calculate left and right trims



83
84
85
86
87
88
89
90
91
92
93
94
 
# File 'lib/genfrag/app/search_command/trim.rb', line 83

def calculate_left_and_right_trims(trim)
  left = {}
  # Should we "dot out" (nucleotide padding) from the primary strand? If no, then we assume the complement needs padding.
  left[:dot_out_from_primary] = (trim[:from_left_primary] > trim[:from_left_complement])
  # How much gets cut off on both primary and complement strands
  left[:trim_from_both] = [trim[:from_left_primary], trim[:from_left_complement]].min

  right = {}
  right[:dot_out_from_primary] = (trim[:from_right_primary] > trim[:from_right_complement])
  right[:trim_from_both] = [trim[:from_right_primary], trim[:from_right_complement]].min
  return [left,right]
end

#calculate_trim_for_nucleotides(re5_ds, re3_ds) ⇒ Object

Keep track of extraneous nucleotides that should be removed from the final fragment

Example BstYI used as RE5 BstYI -

5' - r^g a t c y - 3'
3' - y c t a g^r - 5'

re5_ds.cut_locations.primary              # => [0]
re5_ds.cut_locations.complement           # => [4]
re5_ds.aligned_strands.primary.size       # => 6

# number of nucleotides to trim from the left side on the primary strand
re5_ds.cut_locations.primary.max + 1      # => 1

# number of nucleotides to trim from the left side on the complement strand
re5_ds.cut_locations.complement.max + 1   # => 5

Example BstYI used as RE3 BstYI -

5' - r^g a t c y - 3'
3' - y c t a g^r - 5'

re3_ds.cut_locations.primary              # => [0]
re3_ds.cut_locations.complement           # => [4]
re3_ds.aligned_strands.primary.size       # => 6

# number of nucleotides to trim from the right side on the primary strand
re3_ds.aligned_strands.primary.size - (re3_ds.cut_locations.primary.min + 1)      # => 5

# number of nucleotides to trim from the right side on the complement strand
re3_ds.aligned_strands.primary.size - (re3_ds.cut_locations.complement.min + 1)   # => 1

Example PpiI used as RE5 PpiI -

5' - n n n n n n^n n n n n n n g a a c n n n n n c t c n n n n n n n n n n n n n^n - 3'
3' - n^n n n n n n n n n n n n c t t g n n n n n g a g n n n n n n n n^n n n n n n - 5'

re5_ds.cut_locations.primary              # => [5, 37]
re5_ds.cut_locations.complement           # => [0, 32]
re5_ds.aligned_strands.primary.size       # => 39

# number of nucleotides to trim from the left side on the primary strand
re5_ds.cut_locations.primary.max + 1      # => 38

# number of nucleotides to trim from the left side on the complement strand
re5_ds.cut_locations.complement.max + 1   # => 33

Example PpiI used as RE3 PpiI -

5' - n n n n n n^n n n n n n n g a a c n n n n n c t c n n n n n n n n n n n n n^n - 3'
3' - n^n n n n n n n n n n n n c t t g n n n n n g a g n n n n n n n n^n n n n n n - 5'

re3_ds.cut_locations.primary              # => [5, 37]
re3_ds.cut_locations.complement           # => [0, 32]
re3_ds.aligned_strands.primary.size       # => 39

# number of nucleotides to trim from the right side on the primary strand
re3_ds.aligned_strands.primary.size - (re3_ds.cut_locations.primary.min + 1)      # => 33

# number of nucleotides to trim from the right side on the complement strand
re3_ds.aligned_strands.primary.size - (re3_ds.cut_locations.complement.min + 1)   # => 38



72
73
74
75
76
77
78
79
 
# File 'lib/genfrag/app/search_command/trim.rb', line 72

def calculate_trim_for_nucleotides(re5_ds, re3_ds)
  trim = {}
  trim[:from_left_primary]     = re5_ds.cut_locations.primary.max + 1
  trim[:from_left_complement]  = re5_ds.cut_locations.complement.max + 1
  trim[:from_right_primary]    = re3_ds.aligned_strands.primary.size - (re3_ds.cut_locations.primary.min + 1)
  trim[:from_right_complement] = re3_ds.aligned_strands.primary.size - (re3_ds.cut_locations.complement.min + 1)
  return trim
end

#cli_run(args) ⇒ Object 



6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 
# File 'lib/genfrag/app/search_command.rb', line 6

def cli_run( args )
  parse args

  @input_filenames = ARGV
  input_filenames = [@input_filenames].flatten
  processed_adapters=nil
  
  validate_options(options)
  
  
  if options[:sqlite]
    processed_fasta_file  = SearchCommand::ProcessFile.process_db_fasta_file( SQLite3::Database.new( Genfrag.name_normalized_fasta(input_filenames,options[:filefasta]) + '.db' ) )
    processed_freq_lookup = SearchCommand::ProcessFile.process_db_freq_lookup( SQLite3::Database.new( Genfrag.name_freq_lookup(input_filenames,options[:filefasta],options[:filelookup],options[:re5],options[:re3]) + '.db' ) )
  else
    processed_fasta_file  = SearchCommand::ProcessFile.process_tdf_fasta_file(  IO.readlines( Genfrag.name_normalized_fasta(input_filenames,options[:filefasta]) + '.tdf' ) )
    processed_freq_lookup = SearchCommand::ProcessFile.process_tdf_freq_lookup( IO.readlines( Genfrag.name_freq_lookup(input_filenames,options[:filefasta],options[:filelookup],options[:re5],options[:re3]) + '.tdf' ) )
  end
  
    if options[:fileadapters]
      processed_adapters = SearchCommand::ProcessFile.process_tdf_adapters( IO.readlines( Genfrag.name_adapters(options[:fileadapters]) + '.tdf' ), options[:named_adapter5], options[:named_adapter3] )
    end
  
  run(options, processed_fasta_file, processed_freq_lookup, processed_adapters, true)
end

#left_tail_of(s) ⇒ Object 



120
121
122
123
124
125
126
127
128
129
130
131
132
 
# File 'lib/genfrag/app/search_command/match.rb', line 120

def left_tail_of(s)
# 'PpiI' => "n n n n n n^n n n n n n n g a a c n n n n n c t c n n n n n n n n n n n n n^n"
# => 'nnnnnn'
# 'BstYI' => "r^g a t c y"
# => 'r'

  if s =~ /([^\^]*)\^/
    return $1.tr(' ', '')
  else
    raise "Sequence #{s} has no cuts (defined by symbol '^')"
  end

end

#matches_adapter(five_or_three, primary_frag, complement_frag, raw_frag, trim, primary_strand) ⇒ Object

Does the sequence match the adapter



9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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
 
# File 'lib/genfrag/app/search_command/match.rb', line 9

def matches_adapter(five_or_three, primary_frag, complement_frag, raw_frag, trim, primary_strand)
  adapter_specificity = nil
  adapter_sequence    = nil
  adapter_size        = nil
  trim_primary        = nil
  trim_complement     = nil

  if five_or_three == 5
    tail = right_tail_of(primary_strand)

    adapter_specificity = @adapters[:adapter5_specificity].upcase
    adapter_sequence    = @adapters[:adapter5_sequence].upcase if @adapters[:adapter5_sequence]
    adapter_size        = @adapters[:adapter5_size]
    trim_primary        = trim[:from_left_primary]
    trim_complement     = trim[:from_left_complement]

    # TEMP Check for match
    primary_frag =~ /(\.*)/
    dots_on_primary = $1.size
    lead_in = tail.size + dots_on_primary

    return false if primary_frag[ lead_in .. -1 ].tr('.', '') !~ /^#{adapter_specificity}/i

  elsif five_or_three == 3
    tail = left_tail_of(primary_strand)

    if @adapters[:adapter3_specificity][0].chr == '_'
      adapter_specificity = @adapters[:adapter3_specificity][1..-1].reverse.upcase
    else
      adapter_specificity = Bio::Sequence::NA.new(@adapters[:adapter3_specificity]).forward_complement.to_s.upcase
    end
    adapter_sequence    = Bio::Sequence::NA.new(@adapters[:adapter3_sequence]).forward_complement.to_s.upcase if @adapters[:adapter3_sequence]
    adapter_size        = @adapters[:adapter3_size]
    trim_primary        = trim[:from_right_primary]
    trim_complement     = trim[:from_right_complement]
    primary_frag.reverse!
    complement_frag.reverse!
    raw_frag.reverse!

    # TEMP Check for match
    primary_frag =~ /(\.*)/
    dots_on_primary = $1.size
    lead_in = tail.size + dots_on_primary
    return false if primary_frag[ lead_in .. -1 ].tr('.', '') !~ /^#{adapter_specificity}/i

  else
    raise "First argument to matches_adapter must be a '5' or a '3'. Received: #{five_or_three.inspect}"
  end

  if adapter_sequence
  # adapter-sequence supplied
    new_primary_frag, new_complement_frag = preserve_or_add(adapter_sequence.size, lead_in, adapter_sequence, primary_frag, complement_frag)
  elsif adapter_size
  # adapter-size supplied
    new_primary_frag, new_complement_frag = preserve_or_add(adapter_size, lead_in, adapter_sequence, primary_frag, complement_frag)
  else
  # only the specificity has been provided
    new_primary_frag = ('.' * dots_on_primary) + ('+' * tail.size) + primary_frag[ lead_in .. -1 ]
    new_complement_frag = complement_frag
  end

  if five_or_three == 3
    return [new_primary_frag.reverse, new_complement_frag.reverse]
  else
    return [new_primary_frag, new_complement_frag]
  end
end

#opt_parser ⇒ Object 



31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
 
# File 'lib/genfrag/app/search_command.rb', line 31

def opt_parser
  std_opts = standard_options

  opts = OptionParser.new
  opts.banner = 'Usage: genfrag search [options]'

  opts.separator ''
  opts.separator "  Search a database of sequence fragments that match the last 5'"
  opts.separator "  fragment cut by two restricting enzymes RE3 and RE5, as created by the"
  opts.separator "  index function. Next, adapters are applied to search a subset of"
  opts.separator "  fragments, as is used in some protocols."

  opts.separator ''
  ary = [:verbose, :quiet, :tracktime, :indir, :outdir, :sqlite, :re5, :re3,
    :filelookup, :filefasta, :fileadapters, :adapter5_sequence, :adapter3_sequence,
    :adapter5_size, :adapter3_size, :named_adapter5, :named_adapter3,
    :adapter5, :adapter3, :seqsize
  ]
  ary.each { |a| opts.on(*std_opts[a]) }
      
  opts.separator ''
  opts.separator '  Common Options:'
  opts.on( '-h', '--help', 'show this message' ) { @out.puts opts; exit }
  
  opts.separator '  Examples:'
  opts.separator '    genfrag search -f example.fasta --re5 BstYI --re3 MseI -v'
  opts.separator '    genfrag search -f example.fasta --re5 BstYI --re3 MseI --adapter5 tt'
  opts.separator '    genfrag search -f example.fasta --re5 BstYI --re3 MseI --add 26 --adapter5 ct --adapter3 aa --size 190,215'
  opts.separator '    genfrag search -f example.fasta --re5 BstYI --re3 MseI --adapter5-size 11 --adapter5 tt --adapter3-size 15 --size 168'
  opts.separator '    genfrag search -f example.fasta --re5 BstYI --re3 MseI --adapter5-sequence GACTGCGTAGTGATC --adapter5 tt --size 168'
  opts.separator '    genfrag search -f example.fasta --re5 BstYI --re3 MseI --adapter5-size 11 --adapter5 ct --adapter3-size 15 --adapter3 aa --size 190,215'
  opts.separator '    genfrag search -f example.fasta --re5 BstYI --re3 MseI --add 26 --named-adapter5 BstYI-T4 --named-adapter3 MseI-21 --size 190,215'
  opts
end

#parse(args) ⇒ Object 



66
67
68
69
70
71
72
73
74
75
76
77
 
# File 'lib/genfrag/app/search_command.rb', line 66

def parse( args )
  opts = opt_parser

  if args.empty?
    @out.puts opts
    exit 1
  end
  
# parse the command line arguments
  opts.parse! args
  
end

#preserve_or_add(size, lead_in, adapter_sequence, primary_frag, complement_frag) ⇒ Object 



134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
 
# File 'lib/genfrag/app/search_command/match.rb', line 134

def preserve_or_add(size, lead_in, adapter_sequence, primary_frag, complement_frag)
  if adapter_sequence.nil? or adapter_sequence.empty?
    adapter_sequence = '?' * size
  end
  
  if lead_in >= size
  # need to preserve dots on primary string
    p = ('=' * (lead_in - size)) + adapter_sequence + primary_frag[ lead_in .. -1 ]
    c = complement_frag
  else
  # need to add dots to beginning of complement string
    p = adapter_sequence + primary_frag[ lead_in .. -1 ]
    c = ('=' * (size - lead_in) ) + complement_frag
  end
  [p,c]
end

#right_tail_of(s) ⇒ Object

# Find the fragments that match the search parameters #

def find_matching_fragments(sizes, left, right)
  hits=[]

  s = (@adapters[:adapter5_size] or 0) + (@adapters[:adapter3_size] or 0)

  if [@ops.seqsize].flatten == [0] or [@ops.seqsize].flatten == [nil] or [@ops.seqsize].flatten == ['0']
    sizes.each do |raw_size, info|
      hits << info
    end

  else
    [@ops.seqsize].flatten.each do |seek_size|
      seek_size = seek_size.to_i
      sizes.each do |raw_size, info|
        frag_size = raw_size - left[:trim_from_both] - right[:trim_from_both]
        if (frag_size >= seek_size - s) and (frag_size <= seek_size + s)
          hits << info
        end
      end
    end
  end

  return hits
end



107
108
109
110
111
112
113
114
115
116
117
118
 
# File 'lib/genfrag/app/search_command/match.rb', line 107

def right_tail_of(s)
# 'PpiI' => "n n n n n n^n n n n n n n g a a c n n n n n c t c n n n n n n n n n n n n n^n"
# => 'n'
# 'BstYI' => "r^g a t c y"
# => 'gatcy'

  if s =~ /.*\^(.*)/
    return $1.tr(' ', '')
  else
    raise "Sequence #{s} has no cuts (defined by symbol '^')"
  end
end

#run(ops = OpenStruct.new, processed_fasta_file = nil, processed_freq_lookup = nil, processed_adapters = nil, cli = false) ⇒ Object 



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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
 
# File 'lib/genfrag/app/search_command.rb', line 110

def run(ops=OpenStruct.new, processed_fasta_file=nil, processed_freq_lookup=nil, processed_adapters=nil, cli=false)
  if ops.kind_of? OpenStruct
    @ops = ops.dup
  elsif ops.kind_of? Hash
    @ops = OpenStruct.new(ops)
  else
    raise ArgumentError
  end
  
# Set defaults
  @ops.verbose        ||= false
  @ops.quiet          ||= false
  @ops.sqlite         ||= false
  @ops.re5            ||= nil
  @ops.re3            ||= nil
  @ops.seqsize        ||= [0]
  @ops.adapter5_size  ||= nil
  @ops.adapter3_size  ||= nil
  @ops.adapter5       ||= nil
  @ops.adapter3       ||= nil
  
  @sizes = processed_freq_lookup
  @sequences = processed_fasta_file
  @adapters = {}
  @re5_ds, @re3_ds = [@ops.re5, @ops.re3].map {|x| Bio::RestrictionEnzyme::DoubleStranded.new(x)}
  @re5_ds_aligned_strands_with_cuts_primary = @re5_ds.aligned_strands_with_cuts.primary
  @re3_ds_aligned_strands_with_cuts_primary = @re3_ds.aligned_strands_with_cuts.primary
  if @ops.verbose
    cli_p(cli, <<-END
RE5: #{@ops.re5}
#{@re5_ds_aligned_strands_with_cuts_primary}
#{@re5_ds.aligned_strands_with_cuts.complement}

RE3: #{@ops.re3}
#{@re5_ds_aligned_strands_with_cuts_primary}
#{@re3_ds.aligned_strands_with_cuts.complement}

adapter5: #{@ops.adapter5}
adapter3: #{@ops.adapter3}
END
)
  end

  if @ops.named_adapter5 and @ops.adapter5
    raise ArgumentError, "Cannot have both 'adapter5' and 'named_adapter5'"
  elsif @ops.named_adapter3 and @ops.adapter3
    raise ArgumentError, "Cannot have both 'adapter3' and 'named_adapter3'"
  end
  
  if !processed_adapters and (@ops.named_adapter5 or @ops.named_adapter3)
    raise ArgumentError, "Must specify --fileadapters when using a named_adapter"
  end
  
  if @ops.adapter5_size and @ops.adapter5_sequence
    raise ArgumentError, '--adapter5-sequence and --adapter5-size both supplied, may only have one'
  end
  if @ops.adapter3_size and @ops.adapter3_sequence
    raise ArgumentError, '--adapter3-sequence and --adapter3-size both supplied, may only have one'
  end
  
  if (@ops.adapter5_sequence or @ops.adapter5_size) and !@ops.adapter5
    raise ArgumentError, '--adapter5 missing in presence of --adapter5-sequence or --adapter5-size'
  end
  if (@ops.adapter3_sequence or @ops.adapter3_size) and !@ops.adapter3
    raise ArgumentError, '--adapter3 missing in presence of --adapter3-sequence or --adapter3-size'
  end
  
  if [@ops.seqsize].flatten == [0] or [@ops.seqsize].flatten == [nil] or [@ops.seqsize].flatten == ['0']
    @ops.seqsize = nil
  else
    h = {:ranges => [], :ints => []}
    @ops.seqsize.flatten.each do |s|
      if s.include?('+')
        a = s.split('+')
        c = a[0].to_i
        r = a[1].to_i
        h[:ranges] << (c-r..c+r)
      else
        h[:ints] << s.to_i
      end
    end
    @ops.seqsize = h
  end
  
  if processed_adapters
    adapter_setup_1(processed_adapters)
  else
    adapter_setup_2
  end
  
# translated adapter 3' if given in reverse orientation - e.g. _tt is 
# translated to aa (reversed) and _tct returns the primary strand
# ending in specific 'tct'
  if @adapters[:adapter3_specificity] =~ /^_/
    seq3 = Bio::Sequence::NA.new(@adapters[:adapter3_specificity][1..-1]).downcase
    @adapters[:adapter3_specificity] = seq3.complement.to_s
  end

  @trim = calculate_trim_for_nucleotides(@re5_ds, @re3_ds)
  
# ------
# Start calculations
#
  left_trim, right_trim = calculate_left_and_right_trims(@trim)
  
  results = []
 
  @sizes.values.each do |hit|
    hit.each do |entry|
      seq = @sequences[entry[:fasta_id]][:sequence]
      raw_frag = seq[entry[:offset]..(entry[:offset]+entry[:raw_size]-1)]

      primary_frag, complement_frag = trim_sequences(raw_frag, Bio::Sequence::NA.new(raw_frag).forward_complement, left_trim, right_trim, @trim)
      
      p = primary_frag.dup
      c = complement_frag.dup
      
    # note the next two if-statements at this level chain together with 'p' and 'c'
      if @adapters[:adapter5_specificity]
        p, c = matches_adapter(5, p, c, raw_frag, @trim, @re5_ds_aligned_strands_with_cuts_primary)
        next if !p  # next if returned false -- no match
      end
      
      if @adapters[:adapter3_specificity]
        p, c = matches_adapter(3, p, c, raw_frag, @trim, @re3_ds_aligned_strands_with_cuts_primary)
        next if !p  # next if returned false -- no match
      end
      
      primary_frag_with_adapters = p
      complement_frag_with_adapters = c
              
      if @ops.seqsize
        primary_frag_with_adapters_size = primary_frag_with_adapters.size
        good = false
        if @ops.seqsize[:ints].include?(primary_frag_with_adapters_size)
          good = true
        else
          @ops.seqsize[:ranges].each do |range|
            good = true if range.include?(primary_frag_with_adapters_size)
            break if good
          end
        end
      # next if fragment size not in range
        next if !good
      end
      
      results << {:raw_frag => raw_frag, :primary_frag => primary_frag, :primary_frag_with_adapters => primary_frag_with_adapters, :complement_frag => complement_frag, :complement_frag_with_adapters => complement_frag_with_adapters, :entry => entry, :seq => seq, :definitions => @sequences[entry[:fasta_id]][:definitions]}
    end
  end

  if results.size == 0
    cli_p(cli,"Nothing found") if @ops.verbose
  end

  sorted_results = {}
  results.sort {|a,b| a[:seq] <=> b[:seq]}.each do |r|
    raise "shouldn't happen" if sorted_results[r[:seq]] != nil
    sorted_results[r[:seq]] = {}
    x = sorted_results[r[:seq]]
    x['fasta definition'] = r[:definitions]
    x['sequence size'] = r[:seq].size
    x['fragment - primary strand'] = r[:primary_frag]
    x['fragment - complement strand'] = r[:complement_frag]
    x['fragment with adapters - primary strand'] = r[:primary_frag_with_adapters]
    x['fragment with adapters - complement strand'] = r[:complement_frag_with_adapters]
  end

  if @ops.verbose
    ary = ['fasta definition', 'sequence size', 'fragment - primary strand', 'fragment - complement strand',
      'fragment with adapters - primary strand', 'fragment with adapters - complement strand']    
  else
    ary = ['fragment with adapters - primary strand', 'fragment with adapters - complement strand']    
  end
  sorted_results.each do |k,v|
    cli_p(cli, '---')
    if @ops.verbose
      cli_p(cli, '- sequence')
      cli_p(cli, "  #{k}")
    end

    ary.each do |a|
      cli_p(cli, "- #{a}")
      cli_p(cli, "  #{v[a]}")
    end
  end

  return results
end

#trim_sequences(primary_frag, complement_frag, left, right, trim) ⇒ Object

Do the trimming



98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
 
# File 'lib/genfrag/app/search_command/trim.rb', line 98

def trim_sequences(primary_frag, complement_frag, left, right, trim)
  if left[:dot_out_from_primary]
    primary_frag = "." * trim[:from_left_primary] + primary_frag[trim[:from_left_primary]..-1]
  else
    complement_frag = "." * trim[:from_left_complement] + complement_frag[trim[:from_left_complement]..-1]
  end

  if right[:dot_out_from_primary]
    primary_frag = primary_frag[0..(-1 - trim[:from_right_primary])] + "." * trim[:from_right_primary]
  else
    complement_frag = complement_frag[0..(-1 - trim[:from_right_primary])] + "." * trim[:from_right_primary]
  end

  primary_frag = primary_frag[left[:trim_from_both]..(-1-right[:trim_from_both])]
  complement_frag = complement_frag[left[:trim_from_both]..(-1-right[:trim_from_both])]

  return [primary_frag, complement_frag]
end

#validate_options(o) ⇒ Object 



79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
 
# File 'lib/genfrag/app/search_command.rb', line 79

def validate_options(o)
  if o[:filefasta] == nil
    clierr_p "missing option: must supply fasta filename"
    exit 1
  end
  
  if o[:re5] == nil
    clierr_p "missing option: re5"
    exit 1
  end
  
  if o[:re3] == nil
    clierr_p "missing option: re3"
    exit 1
  end
  
  begin
    Bio::RestrictionEnzyme::DoubleStranded.new(o[:re3])
  rescue
    clierr_p "re3 is not an enzyme name"
    exit 1
  end
  
  begin
    Bio::RestrictionEnzyme::DoubleStranded.new(o[:re5])
  rescue
    clierr_p "re5 is not an enzyme name"
    exit 1
  end
end