Class: ROCker

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

Constant Summary collapse

@@EBIREST =
[ Class ]
'http://www.ebi.ac.uk/Tools'
@@DEFAULTS = 
{
   # General
   :q=>false, :r=>'R', :nucl=>false, :debug=>false,
   # Build
   :positive=>[], :negative=>[], :thr=>2,:genomefrx=>1.0,
	 # ext. software
	 :grinder=>'grinder', :muscle=>'muscle', :blastbins=>'', :seqdepth=>3, :minovl=>0.75,
	 :grindercmd=>'%1$s -reference_file "%2$s" -cf "%3$f" -base_name "%4$s" -dc \'-~*Nn\' -md "uniform 0.1" -mr "95 5" -rd "100 uniform 5"',
	 :musclecmd=>'%1$s -in "%2$s" -out "%3$s" -quiet',
	 :blastcmd=>'%1$s%2$s -query "%3$s" -db "%4$s" -out "%5$s" -num_threads %6$d -outfmt 6 -max_target_seqs 1',
	 :makedbcmd=>'%1$smakeblastdb -dbtype %2$s -in "%3$s" -out "%4$s"',
   # Compile
   :refine=>true, :win=>20, :minscore=>0,
   # Filter
   :sbj=>[],
   # Plot
   :color=>false, :gformat=>'pdf', :width=>9, :height=>9, :impact=>false, :transparency=>true,
}
@@HAS_BUILD_GEMS = 
nil

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(opts) ⇒ ROCker

Returns a new instance of ROCker.



50
51
52
53
54
 
# File 'lib/rocker.rb', line 50

def initialize(opts)
   @o = ROCker.defaults
   opts.each{ |k,v| @o[k] = v }
   RInterface.R_BIN = opts[:r] unless opts[:r].nil?
end

Instance Attribute Details

#o ⇒ Object (readonly)

[ Instance ]


49
50
51
 
# File 'lib/rocker.rb', line 49

def o
  @o
end

Class Method Details

.default(k) ⇒ Object 



35
 
# File 'lib/rocker.rb', line 35

def self.default(k) @@DEFAULTS[k] ; end

.defaults ⇒ Object 



34
 
# File 'lib/rocker.rb', line 34

def self.defaults() @@DEFAULTS ; end

.ebirest ⇒ Object 



33
 
# File 'lib/rocker.rb', line 33

def self.ebirest() @@EBIREST ; end

.has_build_gems? ⇒ Boolean

Returns:

  • (Boolean) 


36
37
38
39
40
41
42
43
44
45
46
 
# File 'lib/rocker.rb', line 36

def self.has_build_gems?
   return @@HAS_BUILD_GEMS unless @@HAS_BUILD_GEMS.nil?
   @@HAS_BUILD_GEMS = TRUE
   begin
	 require 'rubygems'
	 require 'restclient'
   rescue LoadError
	 @@HAS_BUILD_GEMS = FALSE
   end
   @@HAS_BUILD_GEMS
end

Instance Method Details

#bash(cmd, err_msg = nil) ⇒ Object 



468
469
470
471
472
 
# File 'lib/rocker.rb', line 468

def bash(cmd, err_msg=nil)
   o = `#{cmd} 2>&1 && echo '{'`
   raise (err_msg.nil? ? "Error executing: #{cmd}\n\n#{o}" : err_msg) unless o[-2]=='{'
   true
end

#blast2table(blast_f, table_f, aln, minscore) ⇒ Object

[ Utilities ]


425
426
427
428
429
430
431
432
433
434
 
# File 'lib/rocker.rb', line 425

def blast2table(blast_f, table_f, aln, minscore)
   ifh = File.open(blast_f, "r")
   ofh = File.open(table_f, "w")
   while ln = ifh.gets
	 bh = BlastHit.new(ln, aln)
	 ofh.print bh.to_s if bh.bits >= minscore
   end
   ifh.close
   ofh.close
end

#build! ⇒ Object

[ Build ]


57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
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
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
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
 
# File 'lib/rocker.rb', line 57

def build!
   # Check requirements
   puts "Testing environment." unless @o[:q]
   @o[:noblast]=true if @o[:nomg]
   raise "Unsatisfied requirements, please see the help message (-h)." unless ROCker.has_build_gems?
   @o[:positive] += @o[:posori] unless @o[:posori].nil?
   @o[:positive] += File.readlines(@o[:posfile]).map{ |l| l.chomp } unless @o[:posfile].nil?
   @o[:negative] += File.readlines(@o[:negfile]).map{ |l| l.chomp } unless @o[:negfile].nil?
   unless @o[:aln].nil?
      aln = Alignment.new
	 aln.read_fasta @o[:aln]
	 @o[:positive] += aln.get_ids
   end
   raise "-p or -P are mandatory." if @o[:positive].size==0
   raise "-o/--baseout is mandatory." if @o[:baseout].nil?
   if @o[:positive].size == 1 and not @o[:noaln]
	 warn "\nWARNING: Positive set contains only one sequence, turning off alignment.\n\n"
	 @o[:noaln] = true
   end
   self.bash "#{@o[:grinder]} --version", "-G/--grinder must be executable. Is Grinder installed?" unless @o[:nomg]
   self.bash "#{@o[:muscle]} -version", "-M/--muscle must be executable. Is Muscle installed?" unless @o[:noaln]
   self.bash "#{@o[:blastbins]}makeblastdb -version", "-B/--blastbins must contain executables. Is BLAST+ installed?" unless @o[:noblast]
   # Download genes
   puts "Downloading gene data." unless @o[:q]
   f = File.open(@o[:baseout] + '.ref.fasta', 'w')
   if @o[:posori].nil? and @o[:posfile].nil? and not @o[:aln].nil?
	 puts "  * re-using aligned sequences as positive set." unless @o[:q]
	 f.print aln.to_seq_s
	 @o[:noaln] = true
   else
	 puts "  * downloading #{@o[:positive].size} sequence(s) in positive set." unless @o[:q]
	 $stderr.puts "   # #{@o[:positive]}" if @o[:debug]
	 ids = Array.new(@o[:positive])
	 while ids.size>0
  f.print ebiFetch(:uniprotkb, ids.shift(200), :fasta)
	 end
   end
   f.close
   genome_ids = {:positive=>[], :negative=>[]}
   [:positive, :negative].each do |set|
      unless @o[set].size==0
  puts "  * gathering genomes from #{@o[set].size} #{set.to_s} sequence(s)." unless @o[:q]
  $stderr.puts "   # #{@o[set]}" if @o[:debug]
  genome_ids[set] = genes2genomes(@o[set])
	 end
   end
   raise "No genomes associated with the positive set." if genome_ids[:positive].size==0
   genome_ids[:positive] = genome_ids[:positive].sample( (genome_ids[:positive].size*@o[:genomefrx]).round ) if @o[:genomefrx]
   raise "No positive genomes selected for metagenome construction, is --genome-frx too small?" if genome_ids[:positive].empty?
   all_genome_ids = genome_ids.values.reduce(:+).uniq
   
   # Locate genes
   puts "Analyzing genome data." unless @o[:q]
   puts "  * downloading and parsing #{genome_ids[:positive].size} GFF3 document(s)." unless @o[:q]
   $stderr.puts "   # #{genome_ids[:positive]}" if @o[:debug]
   positive_coords = {}
   genome_org = {}
   i = 0
   genome_ids[:positive].each do |genome_id|
	 print "  * scanning #{(i+=1).ordinalize} genome out of #{genome_ids[:positive].size}. \r" unless @o[:q]
	 unless @o[:pertaxon].nil?
  genome_taxon = genome2taxon(genome_id, @o[:pertaxon])
  next unless genome_org[ genome_taxon ].nil?
  genome_org[ genome_taxon ] = genome_id
	 end
	 $stderr.puts "   # Looking for any of #{@o[:positive]}" if @o[:debug]
	 genome_file = @o[:baseout] + '.src.' + i.to_s + '.gff3'
	 if @o[:reuse] and File.exist? genome_file
  puts "  * reusing existing file: #{genome_file}." unless @o[:q]
  ifh = File.open(genome_file, 'r')
  doc = ifh.readlines.grep(/^[^#]/)
  ifh.close
	 else
  genome_file=nil unless @o[:noclean]
  res = ebiFetch(:embl, [genome_id], :gff3, genome_file)
  doc = res.split("\n").grep(/^[^#]/)
	 end
	 doc.each do |ln|
  next if ln =~ /^#/
  r = ln.chomp.split /\t/
  next if r.size < 9
  prots = r[8].split(/;/).grep(/^db_xref=UniProtKB[\/A-Za-z-]*:/){ |xref| xref.split(/:/)[1] }
  p = prots.select{ |p| @o[:positive].include? p }.first
  next if p.nil?
  positive_coords[ r[0] ] ||= []
  positive_coords[ r[0] ] << {
     #:strand	=> r[6],
     :prot_id	=> p,
     :from	=> r[3].to_i,
     :to	=> r[4].to_i
  }
	 end
   end
   print "\n" unless @o[:q]
   unless @o[:pertaxon].nil?
	 genome_ids[:positive] = genome_org.values
	 puts "  Using #{genome_org.size} genome(s) after filtering by #{@o[:pertaxon]}." unless @o[:q]
   end
   all_genome_ids = genome_ids.values.reduce(:+).uniq
   found = positive_coords.values.map{ |a| a.map{ |b| b[:prot_id] } }.reduce(:+)
   raise "Cannot find the genomic location of any provided sequence." if found.nil?
   missing = @o[:positive] - found
   warn "\nWARNING: Cannot find genomic location of sequence(s) #{missing.join(',')}.\n\n" unless missing.size==0 or @o[:genomefrx]<1.0 or not @o[:pertaxon].nil?
   
   # Download genomes
   genomes_file = @o[:baseout] + '.src.fasta'
   if @o[:reuse] and File.exist? genomes_file
	 puts "  * reusing existing file: #{genomes_file}." unless @o[:q]
   else
	 puts "  * downloading #{all_genome_ids.size} genome(s) in FastA." unless @o[:q]
	 $stderr.puts "   # #{all_genome_ids}" if @o[:debug]
	 ids = Array.new(all_genome_ids)
	 ofh = File.open(genomes_file, 'w')
	 while ids.size>0
  ofh.print ebiFetch('embl', ids.shift(200), 'fasta')
	 end
	 ofh.close
   end
   
   # Generate metagenome
   unless @o[:nomg]
	 puts "Generating in silico metagenome" unless @o[:q]
	 if @o[:reuse] and File.exist? @o[:baseout] + ".mg.fasta"
  puts "  * reusing existing file: #{@o[:baseout]}.mg.fasta." unless @o[:q]
	 else
  all_src = File.readlines("#{@o[:baseout]}.src.fasta").select{ |l| l =~ /^>/ }.size
  thrs = [@o[:thr], all_src].min
  puts "  * running grinder and tagging positive reads (#{thrs} threads)." unless @o[:q]
  $stderr.puts "   # #{positive_coords}" if @o[:debug]
  thr_obj = []
  seqs_per_thr = (all_src/thrs).ceil
  (0 .. (thrs-1)).each do |thr_i|
     thr_obj << Thread.new do
 Thread.current[:seqs_a] = thr_i*seqs_per_thr + 1
 Thread.current[:seqs_b] = [Thread.current[:seqs_a] + seqs_per_thr, all_src].min
 # Create sub-fasta
 Thread.current[:ofh] = File.open("#{@o[:baseout]}.src.fasta.#{thr_i.to_s}", 'w')
 Thread.current[:ifh] = File.open("#{@o[:baseout]}.src.fasta", 'r')
 Thread.current[:seq_i] = 0
 while Thread.current[:l] = Thread.current[:ifh].gets
    Thread.current[:seq_i]+=1 if Thread.current[:l] =~ /^>/
    break if Thread.current[:seq_i] > Thread.current[:seqs_b]
    Thread.current[:ofh].print Thread.current[:l] if Thread.current[:seq_i] >= Thread.current[:seqs_a]
 end
 Thread.current[:ifh].close
 Thread.current[:ofh].close
 bash sprintf(@o[:grindercmd], @o[:grinder], "#{@o[:baseout]}.src.fasta.#{thr_i.to_s}", @o[:seqdepth], "#{@o[:baseout]}.mg.tmp.#{thr_i.to_s}")
 # Tag positives
 puts "  * tagging positive reads." unless @o[:q]
 Thread.current[:ifh] = File.open(@o[:baseout] + ".mg.tmp.#{thr_i.to_s}-reads.fa", 'r')
 Thread.current[:ofh] = File.open(@o[:baseout] + ".mg.fasta.#{thr_i.to_s}", 'w')
 while Thread.current[:l]=Thread.current[:ifh].gets
    Thread.current[:rd] = /^>(?<id>\d+) reference=[A-Za-z]+\|(?<genome_id>[A-Za-z0-9_]+)\|.* position=(?<comp>complement\()?(?<from>\d+)\.\.(?<to>\d+)\)? /.match(Thread.current[:l])
    unless Thread.current[:rd].nil?
Thread.current[:positive] = false
positive_coords[Thread.current[:rd][:genome_id]] ||= []
positive_coords[Thread.current[:rd][:genome_id]].each do |gn|
   Thread.current[:left]  = Thread.current[:rd][:to].to_i - gn[:from]
   Thread.current[:right] = gn[:to] - Thread.current[:rd][:from].to_i
   if (Thread.current[:left]*Thread.current[:right] >= 0) and ([Thread.current[:left], Thread.current[:right]].min/(Thread.current[:rd][:to].to_i-Thread.current[:rd][:from].to_i) >= @o[:minovl])
      Thread.current[:positive] = true
      break
   end
end
Thread.current[:l] = ">#{Thread.current[:rd][:id]}#{Thread.current[:positive] ? "@%" : ""} ref=#{Thread.current[:rd][:genome_id]}:#{Thread.current[:rd][:from]}..#{Thread.current[:rd][:to]}#{(Thread.current[:rd][:comp]=='complement(')?'-':'+'}\n"
    end
    Thread.current[:ofh].print Thread.current[:l]
 end
 Thread.current[:ofh].close
 Thread.current[:ifh].close
 Thread.current[:output] = @o[:baseout] + ".mg.fasta.#{thr_i.to_s}"
     end # Thread.new do
  end # (1 .. thrs).each
  # Concatenate results
  ofh = File.open(@o[:baseout] + ".mg.fasta", 'w')
  thr_obj.each do |t|
     t.join
     raise "Thread failed without error trace: #{t}" if t[:output].nil?
     ifh = File.open(t[:output], 'r')
     while l = ifh.gets
        ofh.print l
     end
     ifh.close
     File.unlink t[:output]
  end
  ofh.close
      end
   end # unless @o[:nomg]
   # Align references
   unless @o[:noaln]
	 puts "Aligning reference set." unless @o[:q]
	 if @o[:reuse] and File.exist? "#{@o[:baseout]}.ref.aln"
  puts "  * reusing existing file: #{@o[:baseout]}.ref.aln." unless @o[:q]
	 else
  bash sprintf(@o[:musclecmd], @o[:muscle], "#{@o[:baseout]}.ref.fasta", "#{@o[:baseout]}.ref.aln")
  puts "  +--\n  | IMPORTANT NOTE: Manually checking the alignment before\n  | the 'compile' step is *strongly* encouraged.\n  +--\n" unless @o[:q]
	 end
   end
   # Run BLAST 
   unless @o[:noblast]
	 puts "Running homology search." unless @o[:q]
	 if @o[:reuse] and File.exist? "#{@o[:baseout]}.ref.blast"
  puts "  * reusing existing file: #{@o[:baseout]}.ref.blast." unless @o[:q]
	 else
  puts "  * preparing database." unless @o[:q]
  bash sprintf(@o[:makedbcmd], @o[:blastbins], (@o[:nucl]?'nucl':'prot'), "#{@o[:baseout]}.ref.fasta", "#{@o[:baseout]}.ref")
  puts "  * running BLAST." unless @o[:q]
  bash sprintf(@o[:blastcmd], @o[:blastbins], (@o[:nucl]?'blastn':'blastx'), "#{@o[:baseout]}.mg.fasta", "#{@o[:baseout]}.ref", "#{@o[:baseout]}.ref.blast", @o[:thr])
	 end
   end
   # Clean
   unless @o[:noclean]
	 puts "Cleaning." unless @o[:q]
	 sff  = %w{.src.xml .src.fasta}
	 sff += %w{.mg.tmp-reads.fa .mg.tmp-ranks.txt} unless @o[:nomg]
	 sff += %w{.ref.phr .ref.pin .ref.psq} unless @o[:noblast]
	 sff.each { |sf| File.unlink @o[:baseout] + sf if File.exist? @o[:baseout] + sf }
   end
end

#compile! ⇒ Object

[ Compile ]


278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
 
# File 'lib/rocker.rb', line 278

def compile!
   raise "-a/--alignment is mandatory." if @o[:aln].nil?
   raise "-a/--alignment must exist." unless File.exist? @o[:aln]
   if @o[:table].nil?
	 raise "-t/--table is mandatory unless -b is provided." if @o[:blast].nil?
	 @o[:table] = "#{@o[:blast]}.table"
   end
   raise "-b/--blast is mandatory unless -t exists." if @o[:blast].nil? and not File.exist? @o[:table]
   raise "-k/--rocker is mandatory." if @o[:rocker].nil?

   puts "Testing environment." unless @o[:q]
   bash "echo '' | #{@o[:r]} --vanilla", "-r/--path-to-r must be executable. Is R installed?"
   bash "echo \"library('pROC')\" | #{@o[:r]} --vanilla", "Please install the 'pROC' library for R first."

   puts "Reading files." unless @o[:q]
   puts "  * loading alignment: #{@o[:aln]}." unless @o[:q]
   aln = Alignment.new
   aln.read_fasta @o[:aln]
   
   if File.exist? @o[:table]
	 puts "  * reusing existing file: #{@o[:table]}." unless @o[:q]
   else
	 puts "  * generating table: #{@o[:table]}." unless @o[:q]
	 blast2table(@o[:blast], @o[:table], aln, @o[:minscore])
   end

   puts "Analyzing data." unless @o[:q]
   puts "  * computing windows." unless @o[:q]
   data = ROCData.new(@o[:table], aln, @o[:win])
   data.nucl = @o[:nucl]
   if @o[:refine]
	 puts "  * refining windows." unless @o[:q]
	 warn "Insufficient hits to refine results." unless data.refine! @o[:table]
   end
   puts "  * saving ROCker file: #{@o[:rocker]}." unless @o[:q]
   data.save @o[:rocker]
end

#ebiFetch(db, ids, format, outfile = nil) ⇒ Object 



463
464
465
466
467
 
# File 'lib/rocker.rb', line 463

def ebiFetch(db, ids, format, outfile=nil)
   url = "#{ROCker.ebirest}/dbfetch/dbfetch/#{db.to_s}/#{ids.join(",")}/#{format.to_s}"
   $stderr.puts "   # Calling: #{url}" if @o[:debug]
   self.restcall url
end

#filter! ⇒ Object

[ Filter ]


317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
 
# File 'lib/rocker.rb', line 317

def filter!
   raise "-k/--rocker is mandatory." if @o[:rocker].nil?
   raise "-x/--query-blast is mandatory." if @o[:qblast].nil?
   raise "-o/--out-blast is mandatory." if @o[:oblast].nil?
   
   puts "Reading ROCker file." unless @o[:q]
   data = ROCData.new @o[:rocker]

   puts "Filtering BLAST." unless @o[:q]
   ih = File.open(@o[:qblast], 'r')
   oh = File.open(@o[:oblast], 'w')
   while ln = ih.gets
	 bh = BlastHit.new(ln, data.aln)
	 oh.print ln if not(bh.sfrom.nil?) and bh.bits >= data.win_at_col(bh.midpoint).thr
   end
   ih.close
   oh.close
end

#genes2genomes(gene_ids) ⇒ Object 



435
436
437
438
439
440
441
442
443
 
# File 'lib/rocker.rb', line 435

def genes2genomes(gene_ids)
   genomes = []
   ids = Array.new(gene_ids)
   while ids.size>0
	 doc = ebiFetch(:uniprotkb, ids.shift(200), :annot).split("\n")
	 genomes += doc.grep( /^DR\s+EMBL;/ ).map{ |ln| ln.split('; ')[1] }
   end
   genomes.uniq
end

#genome2taxid(genome_id) ⇒ Object 



444
445
446
447
448
 
# File 'lib/rocker.rb', line 444

def genome2taxid(genome_id)
   ln = ebiFetch('embl', [genome_id], 'annot').split(/[\n\r]/).grep(/^FT\s+\/db_xref="taxon:/).first
   return ln if ln.nil?
   ln.sub(/.*"taxon:(\d+)".*/, "\\1")
end

#genome2taxon(genome_id, rank = 'species') ⇒ Object 



449
450
451
452
 
# File 'lib/rocker.rb', line 449

def genome2taxon(genome_id, rank='species')
   xml = ebiFetch('taxonomy', [genome2taxid(genome_id)], 'enataxonomyxml').gsub(/\s*\n\s*/,'')
   xml.scan(/<taxon [^>]+>/).grep(/rank="#{rank}"/).first.sub(/.* taxId="(\d+)".*/,"\\1")
end

#plot! ⇒ Object

[ Plot ]


344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
 
# File 'lib/rocker.rb', line 344

def plot!
   raise "-k/--rocker is mandatory." if o[:rocker].nil?
   if @o[:table].nil?
	 raise "-t/--table is mandatory unless -b is provided." if @o[:blast].nil?
	 @o[:table] = "#{@o[:blast]}.table"
   end
   raise "-b/--blast is mandatory unless -t exists." if @o[:blast].nil? and not File.exist? @o[:table]

   puts "Testing environment." unless @o[:q]
   bash "echo '' | #{@o[:r]} --vanilla", "-r/--path-to-r must be executable. Is R installed?"

   puts "Reading files." unless @o[:q]
   puts "  * loding ROCker file: #{@o[:rocker]}." unless @o[:q]
   data = ROCData.new @o[:rocker]
   if File.exist? @o[:table]
	 puts "  * reusing existing file: #{@o[:table]}." unless @o[:q]
   else
	 puts "  * generating table: #{@o[:table]}." unless @o[:q]
	 blast2table(@o[:blast], @o[:table], data.aln, @o[:minscore])
   end

   puts "Plotting matches." unless @o[:q]
   extra = @o[:gformat]=='pdf' ? "" : ", units='in', res=300"
   @o[:gout] ||= "#{@o[:rocker]}.#{@o[:gformat]}"
   data.rrun "#{@o[:gformat]}('#{@o[:gout]}', #{@o[:width]}, #{@o[:height]}#{extra});"
   data.rrun "layout(c(2,1,3), heights=c(2-1/#{data.aln.size},3,1));"
   some_thr = data.load_table! @o[:table], @o[:sbj], @o[:minscore]
   data.rrun "par(mar=c(0,4,0,0.5)+.1);"
   data.rrun "plot(1, t='n', xlim=c(0.5,#{data.aln.cols}+0.5), ylim=range(x$V4)+c(-0.04,0.04)*diff(range(x$V4)), xlab='', ylab='Bit score', xaxs='i', xaxt='n');"
   data.rrun "noise <- runif(ncol(x),-.2,.2)"
   data.rrun "arrows(x0=x$V2, x1=x$V3, y0=x$V4+noise, col=ifelse(x$V5==1, rgb(0,0,.5,#{@o[:transparency] ? ".2" : "1"}), rgb(.5,0,0,#{@o[:transparency] ? ".2" : "1"})), length=0);"
   data.rrun "points(x$V6, x$V4+noise, col=ifelse(x$V5==1, rgb(0,0,.5,#{@o[:transparency] ? ".5" : "1"}), rgb(.5,0,0,#{@o[:transparency] ? ".5" : "1"})), pch=19, cex=1/4);"

   puts "Plotting windows." unless @o[:q]
   if some_thr
	 data.rrun "arrows(x0=w$V1, x1=w$V2, y0=w$V5, lwd=2, length=0)"
	 data.rrun "arrows(x0=w$V2[-nrow(w)], x1=w$V1[-1], y0=w$V5[-nrow(w)], y1=w$V5[-1], lwd=2, length=0)"
   end
   data.rrun "legend('bottomright',legend=c('Match span','Match mid-point','Reference','Non-reference')," +
	 "lwd=c(1,NA,1,1),pch=c(NA,19,19,19),col=c('black','black','darkblue','darkred'),ncol=4,bty='n')"

   puts "Plotting alignment." unless @o[:q]
   data.rrun "par(mar=c(0,4,0.5,0.5)+0.1);"
   data.rrun "plot(1, t='n', xlim=c(0,#{data.aln.cols}),ylim=c(1,#{data.aln.seqs.size}),xlab='',ylab='Alignment',xaxs='i',xaxt='n',yaxs='i',yaxt='n',bty='n');"
   i = 0
   data.rrun "clr <- rainbow(26, v=1/2, s=3/4);" if @o[:color]
   data.aln.seqs.values.each do |s|
      color = s.aln.split(//).map{|c| c=="-" ? "'grey80'" : (@o[:sbj].include?(s.id) ? "'red'" : (@o[:color] ? "clr[#{c.ord-64}]" : "'black'"))}.join(',')
	 data.rrun "rect((1:#{data.aln.cols-1})-0.5, rep(#{i}, #{data.aln.cols-1}), (1:#{data.aln.cols-1})+0.5, rep(#{i+1}, #{data.aln.cols-1}), col=c(#{color}), border=NA);"
	 i += 1
   end

   puts "Plotting statistics." unless @o[:q]
   data.rrun "par(mar=c(5,4,0,0.5)+.1);"
   data.rrun "plot(1, t='n', xlim=c(0,#{data.aln.cols}),ylim=c(#{@o[:ylim].nil? ? (@o[:impact] ? "-2,.1" : "50,100") : @o[:ylim]}),xlab='Alignment position (amino acids)',ylab='Precision',xaxs='i');"
   if some_thr
	 sn = data.rrun "100*sum(w$tp)/(sum(w$tp)+sum(w$fn))", :float
	 sp = data.rrun "100*sum(w$tn)/(sum(w$fp)+sum(w$tn))", :float
	 ac = data.rrun "100*(sum(w$tp)+sum(w$tn))/(sum(w$p)+sum(w$n))", :float
	 unless @o[:q]
  puts "  * sensitivity: #{sn}%"
  puts "  * specificity: #{sp}%"
  puts "  * accuracy: #{ac}%"
	 end
	 data.rrun "pos <- (w$V1+w$V2)/2"
	 if @o[:impact]
  data.rrun "lines(pos[!is.na(w$specificity)], (w$specificity[!is.na(w$specificity)]-#{sp})*w$tp[!is.na(w$specificity)]/sum(w$tp), col='darkred', lwd=2, t='o', cex=1/3, pch=19);"
  data.rrun "lines(pos[!is.na(w$sensitivity)], (w$sensitivity[!is.na(w$sensitivity)]-#{sn})*w$tn[!is.na(w$sensitivity)]/sum(w$tn), col='darkgreen', lwd=2, t='o', cex=1/3, pch=19);"
  data.rrun "lines(pos[!is.na(w$accuracy)], (w$accuracy[!is.na(w$accuracy)]-#{ac})*(w$tp+w$tn)[!is.na(w$accuracy)]/sum(c(w$tp, w$tn)), col='darkblue', lwd=2, t='o', cex=1/3, pch=19);"
	 else
  data.rrun "lines(pos[!is.na(w$specificity)], w$specificity[!is.na(w$specificity)], col='darkred', lwd=2, t='o', cex=1/3, pch=19);"
  data.rrun "lines(pos[!is.na(w$sensitivity)], w$sensitivity[!is.na(w$sensitivity)], col='darkgreen', lwd=2, t='o', cex=1/3, pch=19);"
  data.rrun "lines(pos[!is.na(w$accuracy)], w$accuracy[!is.na(w$accuracy)], col='darkblue', lwd=2, t='o', cex=1/3, pch=19);"
	 end
	 #data.rrun "lines(pos[!is.na(w$precision)], w$precision[!is.na(w$precision)], col='purple', lwd=2, t='o', cex=1/3, pch=19);"
   end
   data.rrun "legend('bottomright',legend=c('Specificity','Sensitivity','Accuracy'),lwd=2,col=c('darkred','darkgreen','darkblue'),ncol=3,bty='n')"
   data.rrun "dev.off();"
end

#restcall(url, outfile = nil) ⇒ Object 



453
454
455
456
457
458
459
460
461
462
 
# File 'lib/rocker.rb', line 453

def restcall(url, outfile=nil)
   response = RestClient.get url
   raise "Unable to reach EBI REST client, error code #{response.code}." unless response.code == 200
   unless outfile.nil?
	 ohf = File.open(outfile, 'w')
	 ohf.print response.to_s
	 ohf.close
   end
   response.to_s
end

#search! ⇒ Object

[ Search ]


336
337
338
339
340
341
 
# File 'lib/rocker.rb', line 336

def search!
   raise "-k/--rocker is mandatory." if @o[:rocker].nil?
   raise "Code Under development..."
   # ToDo
   # [ ... ]
end