Tutorial.rd.html 54 KB
 Lucas Nussbaum committed Jun 15, 2011 1 2 3 4 5 6 7 8 9 10 11 12 13  doc/Tutorial.rd

BioRuby Tutorial

• Copyright (C) 2001-2003 KATAYAMA Toshiaki <k .at. bioruby.org>
•  Lucas Nussbaum committed Sep 25, 2011 14 
• Copyright (C) 2005-2011 Pjotr Prins, Naohisa Goto and others
•  Lucas Nussbaum committed Jun 15, 2011 15 
 Lucas Nussbaum committed Sep 25, 2011 16 17 18 

This document was last modified: 2011/03/24 Current editor: Michael O'Keefe <okeefm (at) rpi (dot) edu>

 Lucas Nussbaum committed Jun 15, 2011 19 20 

Introduction

This is a tutorial for using Bioruby. A basic knowledge of Ruby is required.  Lucas Nussbaum committed Sep 25, 2011 21 If you want to know more about the programming language, we recommend the  Lucas Nussbaum committed Jun 15, 2011 22 latest Ruby book Programming Ruby  Lucas Nussbaum committed Sep 25, 2011 23 by Dave Thomas and Andy Hunt - the first edition can be read online  Lucas Nussbaum committed Jun 15, 2011 24 25 26 27 28 here.

For BioRuby you need to install Ruby and the BioRuby package on your computer

You can check whether Ruby is installed on your computer and what version it has with the

% ruby -v
 Lucas Nussbaum committed Sep 25, 2011 29 

command. You should see something like:

 Lucas Nussbaum committed Jun 15, 2011 30 31 32 33 34 35 
ruby 1.8.7 (2008-08-11 patchlevel 72) [i486-linux]

If you see no such thing you'll have to install Ruby using your installation manager. For more information see the Ruby website.

With Ruby download and install Bioruby using the links on the Bioruby website. The recommended installation is via  Lucas Nussbaum committed Sep 25, 2011 36 RubyGems:

 Lucas Nussbaum committed Jun 15, 2011 37 38 39 40 41 42 43 
gem install bio

A lot of BioRuby's documentation exists in the source code and unit tests. To really dive in you will need the latest source code tree. The embedded rdoc documentation can be viewed online at bioruby's rdoc. But first lets start!

Trying Bioruby

 Lucas Nussbaum committed Sep 25, 2011 44 45 46 47 

Bioruby comes with its own shell. After unpacking the sources run one of the following commands:

bioruby

or, from the source tree

cd bioruby 
Lucas Nussbaum committed Jun 15, 2011  48     49     50     51     52     53     54     55     56     57     58     59     60     61     62     63                                                                                                                                                                                                                                                                                                                                                                                                                          ruby -I lib bin/bioruby

and you should see a prompt

bioruby>

Now test the following:

bioruby> require 'bio' bioruby> seq = Bio::Sequence::NA.new("atgcatgcaaaa") ==> "atgcatgcaaaa"  bioruby> seq.complement ==> "ttttgcatgcat"

See the the Bioruby shell section below for more tweaking. If you have trouble running examples also check the section below on trouble shooting. You can also post a question to the mailing list. BioRuby developers usually try to help.

Working with nucleic / amino acid sequences (Bio::Sequence class)

The Bio::Sequence class allows the usual sequence transformations and translations. In the example below the DNA sequence "atgcatgcaaaa" is  Lucas Nussbaum committed Sep 25, 2011 64 65 converted into the complemental strand and spliced into a subsequence; next, the nucleic acid composition is calculated and the sequence is  Lucas Nussbaum committed Jun 15, 2011 66 translated into the amino acid sequence, the molecular weight  Lucas Nussbaum committed Sep 25, 2011 67 68 calculated, and so on. When translating into amino acid sequences, the frame can be specified and optionally the codon table selected (as  Lucas Nussbaum committed Jun 15, 2011 69 70 71 72 73 74 75 76 defined in codontable.rb).

bioruby> seq = Bio::Sequence::NA.new("atgcatgcaaaa") ==> "atgcatgcaaaa"  # complemental sequence (Bio::Sequence::NA object) bioruby> seq.complement ==> "ttttgcatgcat"  
Lucas Nussbaum committed Sep 25, 2011  77                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   bioruby> seq.subseq(3,8) # gets subsequence of positions 3 to 8 (starting from 1) 
Lucas Nussbaum committed Jun 15, 2011  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                                                                                                                                                                                                                                                ==> "gcatgc" bioruby> seq.gc_percent  ==> 33 bioruby> seq.composition  ==> {"a"=>6, "c"=>2, "g"=>2, "t"=>2} bioruby> seq.translate  ==> "MHAK" bioruby> seq.translate(2)        # translate from frame 2 ==> "CMQ" bioruby> seq.translate(1,11)     # codon table 11 ==> "MHAK" bioruby> seq.translate.codes ==> ["Met", "His", "Ala", "Lys"] bioruby> seq.translate.names ==> ["methionine", "histidine", "alanine", "lysine"] bioruby>  seq.translate.composition ==> {"K"=>1, "A"=>1, "M"=>1, "H"=>1} bioruby> seq.translate.molecular_weight ==> 485.605 bioruby> seq.complement.translate ==> "FCMH"

get a random sequence with the same NA count:

bioruby> counts = {'a'=>seq.count('a'),'c'=>seq.count('c'),'g'=>seq.count('g'),'t'=>seq.count('t')} ==> {"a"=>6, "c"=>2, "g"=>2, "t"=>2} bioruby!> randomseq = Bio::Sequence::NA.randomize(counts)  ==!> "aaacatgaagtc"  bioruby!> print counts a6c2g2t2   bioruby!> p counts {"a"=>6, "c"=>2, "g"=>2, "t"=>2}

The p, print and puts methods are standard Ruby ways of outputting to the screen. If you want to know more about standard Ruby commands you can use the 'ri' command on the command line (or the help command in Windows). For example

% ri puts % ri p % ri File.open
 Lucas Nussbaum committed Sep 25, 2011 116 117 

Nucleic acid sequence are members of the Bio::Sequence::NA class, and amino acid sequence are members of the Bio::Sequence::AA class. Shared  Lucas Nussbaum committed Jun 15, 2011 118 methods are in the parent Bio::Sequence class.

 Lucas Nussbaum committed Sep 25, 2011 119 

As Bio::Sequence inherits Ruby's String class, you can use  Lucas Nussbaum committed Jun 15, 2011 120 121 122 123 124 125 126 127 128 129 130 131 String class methods. For example, to get a subsequence, you can not only use subseq(from, to) but also String#[].

Please take note that the Ruby's string's are base 0 - i.e. the first letter has index 0, for example:

bioruby> s = 'abc' ==> "abc" bioruby> s[0].chr ==> "a" bioruby> s[0..1] ==> "ab"

So when using String methods, you should subtract 1 from positions conventionally used in biology. (subseq method will throw an exception if you  Lucas Nussbaum committed Sep 25, 2011 132 specify positions smaller than or equal to 0 for either one of the "from" or "to".)

 Lucas Nussbaum committed Jun 15, 2011 133 134 135 136 137 

The window_search(window_size, step_size) method shows a typical Ruby way of writing concise and clear code using 'closures'. Each sliding window creates a subsequence which is supplied to the enclosed block through a variable named +s+.

 Lucas Nussbaum committed Sep 25, 2011 138 
• Show average percentage of GC content for 20 bases (stepping the default one base at a time):

 Lucas Nussbaum committed Jun 15, 2011 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 
bioruby> seq = Bio::Sequence::NA.new("atgcatgcaattaagctaatcccaattagatcatcccgatcatcaaaaaaaaaa") ==> "atgcatgcaattaagctaatcccaattagatcatcccgatcatcaaaaaaaaaa"  bioruby> a=[]; seq.window_search(20) { |s| a.push s.gc_percent }  bioruby> a ==> [30, 35, 40, 40, 35, 35, 35, 30, 25, 30, 30, 30, 35, 35, 35, 35, 35, 40, 45, 45, 45, 45, 40, 35, 40, 40, 40, 40, 40, 35, 35, 35, 30, 30, 30]

Since the class of each subsequence is the same as original sequence (Bio::Sequence::NA or Bio::Sequence::AA or Bio::Sequence), you can use all methods on the subsequence. For example,

• Shows translation results for 15 bases shifting a codon at a time

Finally, the window_search method returns the last leftover subsequence. This allows for example

• Divide a genome sequence into sections of 10000bp and output FASTA formatted sequences (line width 60 chars). The 1000bp at the start and end of each subsequence overlapped. At the 3' end of the sequence the leftover is also added:

i = 1 textwidth=60 remainder = seq.window_search(10000, 9000) do |s|   puts s.to_fasta("segment #{i}", textwidth)   i += 1 end if remainder   puts remainder.to_fasta("segment #{i}", textwidth)  end

If you don't want the overlapping window, set window size and stepping size to equal values.

Other examples

• Count the codon usage

bioruby> codon_usage = Hash.new(0) bioruby> seq.window_search(3, 3) { |s| codon_usage[s] += 1 } bioruby> codon_usage ==> {"cat"=>1, "aaa"=>3, "cca"=>1, "att"=>2, "aga"=>1, "atc"=>1, "cta"=>1, "gca"=>1, "cga"=>1, "tca"=>3, "aag"=>1, "tcc"=>1, "atg"=>1}
• Calculate molecular weight for each 10-aa peptide (or 10-nt nucleic acid)

bioruby> a = [] bioruby> seq.window_search(10, 10) { |s| a.push s.molecular_weight } bioruby> a ==> [3096.2062, 3086.1962, 3056.1762, 3023.1262, 3073.2262]

require 'bio'  input_seq = ARGF.read       # reads all files in arguments  my_naseq = Bio::Sequence::NA.new(input_seq) my_aaseq = my_naseq.translate  puts my_aaseq
 Lucas Nussbaum committed Sep 25, 2011 198 199 

Save the program above as na2aa.rb. Prepare a nucleic acid sequence described below and save it as my_naseq.txt:

 Lucas Nussbaum committed Jun 15, 2011 200 201 202 203 204 205 206 207 208 209 
gtggcgatctttccgaaagcgatgactggagcgaagaaccaaagcagtgacatttgtctg atgccgcacgtaggcctgataagacgcggacagcgtcgcatcaggcatcttgtgcaaatg tcggatgcggcgtga

na2aa.rb translates a nucleic acid sequence to a protein sequence. For example, translates my_naseq.txt:

% ruby na2aa.rb my_naseq.txt

or use a pipe!

% cat my_naseq.txt|ruby na2aa.rb

Outputs

VAIFPKAMTGAKNQSSDICLMPHVGLIRRGQRRIRHLVQMSDAA*
 Lucas Nussbaum committed Sep 25, 2011 210 

You can also write this, a bit fancifully, as a one-liner script.

 Lucas Nussbaum committed Jun 15, 2011 211 212 213 214 215 216 217 
% ruby -r bio -e 'p Bio::Sequence::NA.new(<.read).translate' my_naseq.txt In the next section we will retrieve data from databases instead of using raw sequence files. One generic example of the above can be found in ./sample/na2aa.rb. Parsing GenBank data (Bio::GenBank class) We assume that you already have some GenBank data files. (If you don't, download some .seq files from ftp://ftp.ncbi.nih.gov/genbank/)  Lucas Nussbaum committed Sep 25, 2011 218  As an example we will fetch the ID, definition and sequence of each entry  Lucas Nussbaum committed Jun 15, 2011 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 from the GenBank format and convert it to FASTA. This is also an example script in the BioRuby distribution. A first attempt could be to use the Bio::GenBank class for reading in the data: #!/usr/bin/env ruby require 'bio' # Read all lines from STDIN split by the GenBank delimiter while entry = gets(Bio::GenBank::DELIMITER) gb = Bio::GenBank.new(entry) # creates GenBank object print ">#{gb.accession} " # Accession puts gb.definition # Definition puts gb.naseq # Nucleic acid sequence # (Bio::Sequence::NA object) end But that has the disadvantage the code is tied to GenBank input. A more generic method is to use Bio::FlatFile which allows you to use different input formats: #!/usr/bin/env ruby require 'bio' ff = Bio::FlatFile.new(Bio::GenBank, ARGF) ff.each_entry do |gb| definition = "#{gb.accession} #{gb.definition}" puts gb.naseq.to_fasta(definition, 60) end For example, in turn, reading FASTA format files: #!/usr/bin/env ruby require 'bio' ff = Bio::FlatFile.new(Bio::FastaFormat, ARGF) ff.each_entry do |f| puts "definition : " + f.definition puts "nalen : " + f.nalen.to_s puts "naseq : " + f.naseq end  Lucas Nussbaum committed Sep 25, 2011 259  In the above two scripts, the first arguments of Bio::FlatFile.new are  Lucas Nussbaum committed Jun 15, 2011 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 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 database classes of BioRuby. This is expanded on in a later section. Again another option is to use the Bio::DB.open class: #!/usr/bin/env ruby require 'bio' ff = Bio::GenBank.open("gbvrl1.seq") ff.each_entry do |gb| definition = "#{gb.accession} #{gb.definition}" puts gb.naseq.to_fasta(definition, 60) end Next, we are going to parse the GenBank 'features', which is normally very complicated: #!/usr/bin/env ruby require 'bio' ff = Bio::FlatFile.new(Bio::GenBank, ARGF) # iterates over each GenBank entry ff.each_entry do |gb| # shows accession and organism puts "# #{gb.accession} - #{gb.organism}" # iterates over each element in 'features' gb.features.each do |feature| position = feature.position hash = feature.assoc # put into Hash # skips the entry if "/translation=" is not found next unless hash['translation'] # collects gene name and so on and joins it into a string gene_info = [ hash['gene'], hash['product'], hash['note'], hash['function'] ].compact.join(', ') # shows nucleic acid sequence puts ">NA splicing('#{position}') : #{gene_info}" puts gb.naseq.splicing(position) # shows amino acid sequence translated from nucleic acid sequence puts ">AA translated by splicing('#{position}').translate" puts gb.naseq.splicing(position).translate # shows amino acid sequence in the database entry (/translation=) puts ">AA original translation" puts hash['translation'] end end • Note: In this example Feature#assoc method makes a Hash from a feature object. It is useful because you can get data from the hash  Lucas Nussbaum committed Sep 25, 2011 314  by using qualifiers as keys. But there is a risk some information is lost when two or more qualifiers are the same. Therefore an Array is returned by Feature#feature. •  Lucas Nussbaum committed Jun 15, 2011 315   Lucas Nussbaum committed Sep 25, 2011 316  Bio::Sequence#splicing splices subsequences from nucleic acid sequences  Lucas Nussbaum committed Jun 15, 2011 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 according to location information used in GenBank, EMBL and DDBJ. When the specified translation table is different from the default (universal), or when the first codon is not "atg" or the protein contains selenocysteine, the two amino acid sequences will differ. The Bio::Sequence#splicing method takes not only DDBJ/EMBL/GenBank feature style location text but also Bio::Locations object. For more information about location format and Bio::Locations class, see bio/location.rb. • Splice according to location string used in a GenBank entry naseq.splicing('join(2035..2050,complement(1775..1818),13..345') • Generate Bio::Locations object and pass the splicing method locs = Bio::Locations.new('join((8298.8300)..10206,1..855)') naseq.splicing(locs)  Lucas Nussbaum committed Sep 25, 2011 332  You can also use this splicing method for amino acid sequences  Lucas Nussbaum committed Jun 15, 2011 333 334 335 336 337 338 339 340 341 342 343 (Bio::Sequence::AA objects). • Splicing peptide from a protein (e.g. signal peptide) aaseq.splicing('21..119') More databases Databases in BioRuby are essentially accessed like that of GenBank with classes like Bio::GenBank, Bio::KEGG::GENES. A full list can be found in the ./lib/bio/db directory of the BioRuby source tree. In many cases the Bio::DatabaseClass acts as a factory pattern and recognises the database type automatically - returning a  Lucas Nussbaum committed Sep 25, 2011 344 parsed object. For example using Bio::FlatFile class as described above. The first argument of the Bio::FlatFile.new is database class name in BioRuby (such as Bio::GenBank, Bio::KEGG::GENES and so on).  Lucas Nussbaum committed Jun 15, 2011 345 346 347 348 349 350 351 352 353 354 355 356 357  ff = Bio::FlatFile.new(Bio::DatabaseClass, ARGF) Isn't it wonderful that Bio::FlatFile automagically recognizes each database class? #!/usr/bin/env ruby require 'bio' ff = Bio::FlatFile.auto(ARGF) ff.each_entry do |entry| p entry.entry_id # identifier of the entry p entry.definition # definition of the entry p entry.seq # sequence data of the entry end  Lucas Nussbaum committed Sep 25, 2011 358  An example that can take any input, filter using a regular expression and output  Lucas Nussbaum committed Jun 15, 2011 359 360 361 362 to a FASTA file can be found in sample/any2fasta.rb. With this technique it is possible to write a Unix type grep/sort pipe for sequence information. One example using scripts in the BIORUBY sample folder: fastagrep.rb '/At|Dm/' database.seq | fastasort.rb  Lucas Nussbaum committed Sep 25, 2011 363  greps the database for Arabidopsis and Drosophila entries and sorts the output to FASTA.  Lucas Nussbaum committed Jun 15, 2011 364 365  Other methods to extract specific data from database objects can be different between databases, though some methods are common (see the  Lucas Nussbaum committed Sep 25, 2011 366 guidelines for common methods in bio/db.rb).  Lucas Nussbaum committed Jun 15, 2011 367 368 369 370 371 372 373 374 375  • entry_id --> gets ID of the entry • definition --> gets definition of the entry • reference --> gets references as Bio::Reference object • organism --> gets species • seq, naseq, aaseq --> returns sequence as corresponding sequence object Refer to the documents of each database to find the exact naming of the included methods.  Lucas Nussbaum committed Sep 25, 2011 376 377  In general, BioRuby uses the following conventions: when a method name is plural, the method returns some object as an Array. For  Lucas Nussbaum committed Jun 15, 2011 378 379 380 381 example, some classes have a "references" method which returns multiple Bio::Reference objects as an Array. And some classes have a "reference" method which returns a single Bio::Reference object. Alignments (Bio::Alignment)  Lucas Nussbaum committed Sep 25, 2011 382  The Bio::Alignment class in bio/alignment.rb is a container class like Ruby's Hash and Array classes and BioPerl's Bio::SimpleAlign. A very simple example is:  Lucas Nussbaum committed Jun 15, 2011 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  bioruby> seqs = [ 'atgca', 'aagca', 'acgca', 'acgcg' ] bioruby> seqs = seqs.collect{ |x| Bio::Sequence::NA.new(x) } # creates alignment object bioruby> a = Bio::Alignment.new(seqs) bioruby> a.consensus ==> "a?gc?" # shows IUPAC consensus p a.consensus_iupac # ==> "ahgcr" # iterates over each seq a.each { |x| p x } # ==> # "atgca" # "aagca" # "acgca" # "acgcg" # iterates over each site a.each_site { |x| p x } # ==> # ["a", "a", "a", "a"] # ["t", "a", "c", "c"] # ["g", "g", "g", "g"] # ["c", "c", "c", "c"] # ["a", "a", "a", "g"] # doing alignment by using CLUSTAL W. # clustalw command must be installed. factory = Bio::ClustalW.new a2 = a.do_align(factory)  Lucas Nussbaum committed Sep 25, 2011 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428  Read a ClustalW or Muscle 'ALN' alignment file: bioruby> aln = Bio::ClustalW::Report.new(File.read('../test/data/clustalw/example1.aln')) bioruby> aln.header ==> "CLUSTAL 2.0.9 multiple sequence alignment" Fetch a sequence: bioruby> seq = aln.get_sequence(1) bioruby> seq.definition ==> "gi|115023|sp|P10425|" Get a partial sequence: bioruby> seq.to_s[60..120] ==> "LGYFNG-EAVPSNGLVLNTSKGLVLVDSSWDNKLTKELIEMVEKKFQKRVTDVIITHAHAD" Show the full alignment residue match information for the sequences in the set: bioruby> aln.match_line[60..120] ==> " . **. . .. ::*: . * : : . .: .* * *" Return a Bio::Alignment object: bioruby> aln.alignment.consensus[60..120] ==> "???????????SN?????????????D??????????L??????????????????H?H?D"  Lucas Nussbaum committed Jun 15, 2011 429 430 431  Restriction Enzymes (Bio::RE) BioRuby has extensive support for restriction enzymes (REs). It contains a full library of commonly used REs (from REBASE) which can be used to cut single  Lucas Nussbaum committed Sep 25, 2011 432 stranded RNA or double stranded DNA into fragments. To list all enzymes:  Lucas Nussbaum committed Jun 15, 2011 433 434 435 436  rebase = Bio::RestrictionEnzyme.rebase rebase.each do |enzyme_name, info| p enzyme_name end  Lucas Nussbaum committed Sep 25, 2011 437  and to cut a sequence with an enzyme follow up with:  Lucas Nussbaum committed Jun 15, 2011 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462  res = seq.cut_with_enzyme('EcoRII', {:max_permutations => 0}, {:view_ranges => true}) if res.kind_of? Symbol #error err = Err.find_by_code(res.to_s) unless err err = Err.new(:code => res.to_s) end end res.each do |frag| em = EnzymeMatch.new em.p_left = frag.p_left em.p_right = frag.p_right em.c_left = frag.c_left em.c_right = frag.c_right em.err = nil em.enzyme = ar_enz em.sequence = ar_seq p em end Sequence homology search by using the FASTA program (Bio::Fasta) Let's start with a query.pep file which contains a sequence in FASTA format. In this example we are going to execute a homology search from a remote internet site or on your local machine. Note that you  Lucas Nussbaum committed Sep 25, 2011 463 can use the ssearch program instead of fasta when you use it in your  Lucas Nussbaum committed Jun 15, 2011 464 465 466 local machine. using FASTA in local machine Install the fasta program on your machine (the command name looks like  Lucas Nussbaum committed Sep 25, 2011 467 468 fasta34. FASTA can be downloaded from ftp://ftp.virginia.edu/pub/fasta/). First, you must prepare your FASTA-formatted database sequence file  Lucas Nussbaum committed Jun 15, 2011 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 target.pep and FASTA-formatted query.pep. #!/usr/bin/env ruby require 'bio' # Creates FASTA factory object ("ssearch" instead of # "fasta34" can also work) factory = Bio::Fasta.local('fasta34', ARGV.pop) (EDITOR's NOTE: not consistent pop command) ff = Bio::FlatFile.new(Bio::FastaFormat, ARGF) # Iterates over each entry. the variable "entry" is a # Bio::FastaFormat object: ff.each do |entry| # shows definition line (begins with '>') to the standard error outputstderr.puts "Searching ... " + entry.definition    # executes homology search. Returns Bio::Fasta::Report object.   report = factory.query(entry)    # Iterates over each hit   report.each do |hit|     # If E-value is smaller than 0.0001     if hit.evalue < 0.0001       # shows identifier of query and hit, E-value, start and        # end positions of homologous region        print "#{hit.query_id} : evalue #{hit.evalue}\t#{hit.target_id} at "       p hit.lap_at     end   end end
 Lucas Nussbaum committed Sep 25, 2011 501 

We named above script f_search.rb. You can execute it as follows:

 Lucas Nussbaum committed Jun 15, 2011 502 503 504 505 506 507 
% ./f_search.rb query.pep target.pep > f_search.out

In above script, the variable "factory" is a factory object for executing FASTA many times easily. Instead of using Fasta#query method, Bio::Sequence#fasta method can be used.

seq = ">test seq\nYQVLEEIGRGSFGSVRKVIHIPTKKLLVRKDIKYGHMNSKE" seq.fasta(factory)
 Lucas Nussbaum committed Sep 25, 2011 508 509 

When you want to add options to FASTA commands, you can set the third argument of the Bio::Fasta.local method. For example, the following sets ktup to 1 and gets a list of the top 10 hits:

 Lucas Nussbaum committed Jun 15, 2011 510 511 
factory = Bio::Fasta.local('fasta34', 'target.pep', '-b 10') factory.ktup = 1
 Lucas Nussbaum committed Sep 25, 2011 512 

Bio::Fasta#query returns a Bio::Fasta::Report object.  Lucas Nussbaum committed Jun 15, 2011 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 We can get almost all information described in FASTA report text with the Report object. For example, getting information for hits:

report.each do |hit|   puts hit.evalue           # E-value   puts hit.sw               # Smith-Waterman score (*)   puts hit.identity         # % identity   puts hit.overlap          # length of overlapping region   puts hit.query_id         # identifier of query sequence   puts hit.query_def        # definition(comment line) of query sequence   puts hit.query_len        # length of query sequence   puts hit.query_seq        # sequence of homologous region   puts hit.target_id        # identifier of hit sequence   puts hit.target_def       # definition(comment line) of hit sequence   puts hit.target_len       # length of hit sequence   puts hit.target_seq       # hit of homologous region of hit sequence   puts hit.query_start      # start position of homologous                              # region in query sequence   puts hit.query_end        # end position of homologous region                              # in query sequence   puts hit.target_start     # start posiotion of homologous region                              # in hit(target) sequence   puts hit.target_end       # end position of homologous region                              # in hit(target) sequence   puts hit.lap_at           # array of above four numbers end
 Lucas Nussbaum committed Sep 25, 2011 538 539 

Most of above methods are common to the Bio::Blast::Report described below. Please refer to the documentation of the Bio::Fasta::Report class for  Lucas Nussbaum committed Jun 15, 2011 540 FASTA-specific details.

 Lucas Nussbaum committed Sep 25, 2011 541 

If you need the original output text of FASTA program you can use the "output" method of the factory object after the "query" method.

 Lucas Nussbaum committed Jun 15, 2011 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 
report = factory.query(entry) puts factory.output

using FASTA from a remote internet site

• Note: Currently, only GenomeNet (fasta.genome.jp) is supported. check the class documentation for updates.

For accessing a remote site the Bio::Fasta.remote method is used instead of Bio::Fasta.local. When using a remote method, the databases available may be limited, but, otherwise, you can do the same things as with a local method.

Available databases in GenomeNet:

• Protein database
• nr-aa, genes, vgenes.pep, swissprot, swissprot-upd, pir, prf, pdbstr
• Nucleic acid database
• nr-nt, genbank-nonst, gbnonst-upd, dbest, dbgss, htgs, dbsts, embl-nonst, embnonst-upd, genes-nt, genome, vgenes.nuc

Select the databases you require. Next, give the search program from the type of query sequence and database.

 Lucas Nussbaum committed Sep 25, 2011 568 
• When query is an amino acid sequence  Lucas Nussbaum committed Jun 15, 2011 569 570 571 572 573 574 575 
• When protein database, program is "fasta".
• When nucleic database, program is "tfasta".
• When query is a nucleic acid sequence
• When nucleic database, program is "fasta".
•  Lucas Nussbaum committed Sep 25, 2011 576 
• (When protein database, the search would fail.)
•  Lucas Nussbaum committed Jun 15, 2011 577 578 
 Lucas Nussbaum committed Sep 25, 2011 579 

For example, run:

 Lucas Nussbaum committed Jun 15, 2011 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 
program = 'fasta' database = 'genes'  factory = Bio::Fasta.remote(program, database)

and try out the same commands as with the local search shown earlier.

Homology search by using BLAST (Bio::Blast class)

The BLAST interface is very similar to that of FASTA and both local and remote execution are supported. Basically replace above examples Bio::Fasta with Bio::Blast!

For example the BLAST version of f_search.rb is:

# create BLAST factory object factory = Bio::Blast.local('blastp', ARGV.pop)

For remote execution of BLAST in GenomeNet, Bio::Blast.remote is used. The parameter "program" is different from FASTA - as you can expect:

• When query is a amino acid sequence
• When protein database, program is "blastp".
• When nucleic database, program is "tblastn".
• When query is a nucleic acid sequence
• When protein database, program is "blastx"
• When nucleic database, program is "blastn".
• ("tblastx" for six-frame search.)

Bio::BLAST uses "-m 7" XML output of BLAST by default when either XMLParser or REXML (both of them are XML parser libraries for Ruby - of the two XMLParser is the fastest) is installed on your computer. In  Lucas Nussbaum committed Sep 25, 2011 610 Ruby version 1.8.0 or later, REXML is bundled with Ruby's  Lucas Nussbaum committed Jun 15, 2011 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 distribution.

When no XML parser library is present, Bio::BLAST uses "-m 8" tabular deliminated format. Available information is limited with the "-m 8" format so installing an XML parser is recommended.

Again, the methods in Bio::Fasta::Report and Bio::Blast::Report (and Bio::Fasta::Report::Hit and Bio::Blast::Report::Hit) are similar. There are some additional BLAST methods, for example, bit_score and midline.

report.each do |hit|   puts hit.bit_score          puts hit.query_seq          puts hit.midline            puts hit.target_seq          puts hit.evalue             puts hit.identity           puts hit.overlap            puts hit.query_id           puts hit.query_def          puts hit.query_len          puts hit.target_id          puts hit.target_def         puts hit.target_len         puts hit.query_start        puts hit.query_end          puts hit.target_start       puts hit.target_end         puts hit.lap_at           end

For simplicity and API compatibility, some information such as score  Lucas Nussbaum committed Sep 25, 2011 641 is extracted from the first Hsp (High-scoring Segment Pair).

 Lucas Nussbaum committed Jun 15, 2011 642 

Check the documentation for Bio::Blast::Report to see what can be  Lucas Nussbaum committed Sep 25, 2011 643 retrieved. For now suffice to say that Bio::Blast::Report has a  Lucas Nussbaum committed Jun 15, 2011 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 hierarchical structure mirroring the general BLAST output stream:

• In a Bio::Blast::Report object, @iterations is an array of Bio::Blast::Report::Iteration objects.
• In a Bio::Blast::Report::Iteration object, @hits is an array of Bio::Blast::Report::Hits objects.
• In a Bio::Blast::Report::Hits object, @hsps is an array of Bio::Blast::Report::Hsp objects.

Parsing existing BLAST output files

When you already have BLAST output files and you want to parse them, you can directly create Bio::Blast::Report objects without the Bio::Blast factory object. For this purpose use Bio::Blast.reports, which supports the "-m 0" default and "-m 7" XML type output format.

• For example:

blast_version = nil; result = [] Bio::Blast.reports(File.new("../test/data/blast/blastp-multi.m7")) do |report|   blast_version = report.version   report.iterations.each do |itr|     itr.hits.each do |hit|       result.push hit.target_id     end   end end blast_version # ==> "blastp 2.2.18 [Mar-02-2008]" result # ==> ["BAB38768", "BAB38768", "BAB38769", "BAB37741"]
• another example:

require 'bio' Bio::Blast.reports(ARGF) do |report|    puts "Hits for " + report.query_def + " against " + report.db   report.each do |hit|     print hit.target_id, "\t", hit.evalue, "\n" if hit.evalue < 0.001   end end

Save the script as hits_under_0.001.rb and to process BLAST output files *.xml, you can run it with:

% ruby hits_under_0.001.rb *.xml

Sometimes BLAST XML output may be wrong and can not be parsed. Check whether blast is version 2.2.5 or later. See also blast --help.

Bio::Blast loads the full XML file into memory. If this causes a problem you can split the BLAST XML file into smaller chunks using XML-Twig. An example can be found in Biotools.

Note: this section is an advanced topic

Here a more advanced application for using BLAST sequence homology search services. BioRuby currently only supports GenomeNet. If you want to add other sites, you must write the following:

• the calling CGI (command-line options must be processed for the site).
• make sure you get BLAST output text as supported format by BioRuby (e.g. "-m 8", "-m 7" or default("-m 0")).

In addition, you must write a private class method in Bio::Blast named "exec_MYSITE" to get query sequence and to pass the result to Bio::Blast::Report.new(or Bio::Blast::Default::Report.new):

factory = Bio::Blast.remote(program, db, option, 'MYSITE')
 Lucas Nussbaum committed Sep 25, 2011 709 

When you write above routines, please send them to the BioRuby project, and they may be included in future releases.

 Lucas Nussbaum committed Jun 15, 2011 710 711 

Generate a reference list using PubMed (Bio::PubMed)

Nowadays using NCBI E-Utils is recommended. Use Bio::PubMed.esearch  Lucas Nussbaum committed Sep 25, 2011 712 and Bio::PubMed.efetch.

 Lucas Nussbaum committed Jun 15, 2011 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 
#!/usr/bin/env ruby  require 'bio'  # NCBI announces that queries without email address will return error # after June 2010. When you modify the script, please enter your email # address instead of the staff's. Bio::NCBI.default_email = 'staff@bioruby.org'  keywords = ARGV.join(' ')  options = {   'maxdate' => '2003/05/31',   'retmax' => 1000, }  entries = Bio::PubMed.esearch(keywords, options)  Bio::PubMed.efetch(entries).each do |entry|   medline = Bio::MEDLINE.new(entry)   reference = medline.reference   puts reference.bibtex end

The script works same as pmsearch.rb. But, by using NCBI E-Utils, more options are available. For example published dates to search and maximum number of hits to show results can be specified.

See the help page of E-Utils for more details.

In this section, we explain the simple usage of TeX for the BibTeX format bibliography list collected by above scripts. For example, to save BibTeX format bibliography data to a file named genoinfo.bib.

% ./pmfetch.rb 10592173 >> genoinfo.bib % ./pmsearch.rb genome bioinformatics >> genoinfo.bib

The BibTeX can be used with Tex or LaTeX to form bibliography information with your journal article. For more information  Lucas Nussbaum committed Sep 25, 2011 750 on using BibTex see BibTex HowTo site. A quick example:

 Lucas Nussbaum committed Jun 15, 2011 751 752 753 754 755 756 757 758 759 760 761 762 

Save this to hoge.tex:

\documentclass{jarticle} \begin{document} \bibliographystyle{plain} foo bar KEGG database~\cite{PMID:10592173} baz hoge fuga. \bibliography{genoinfo} \end{document}

Then,

% latex hoge % bibtex hoge # processes genoinfo.bib % latex hoge  # creates bibliography list % latex hoge  # inserts correct bibliography reference
 Lucas Nussbaum committed Sep 25, 2011 763 

Now, you get hoge.dvi and hoge.ps - the latter of which can be viewed with any Postscript viewer.

 Lucas Nussbaum committed Jun 15, 2011 764 765 766 

Bio::Reference#bibitem

When you don't want to create a bib file, you can use Bio::Reference#bibitem method instead of Bio::Reference#bibtex.  Lucas Nussbaum committed Sep 25, 2011 767 In the above pmfetch.rb and pmsearch.rb scripts, change

 Lucas Nussbaum committed Jun 15, 2011 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 
puts reference.bibtex

to

puts reference.bibitem

Output documents should be bundled in \begin{thebibliography} and \end{thebibliography}. Save the following to hoge.tex

\documentclass{jarticle} \begin{document} foo bar KEGG database~\cite{PMID:10592173} baz hoge fuga.  \begin{thebibliography}{00}  \bibitem{PMID:10592173} Kanehisa, M., Goto, S. KEGG: kyoto encyclopedia of genes and genomes., {\em Nucleic Acids Res}, 28(1):27--30, 2000.  \end{thebibliography} \end{document}

and run

% latex hoge   # creates bibliography list % latex hoge   # inserts corrent bibliography reference

OBDA

OBDA (Open Bio Database Access) is a standardized method of sequence database access developed by the Open Bioinformatics Foundation. It was created during the BioHackathon by BioPerl, BioJava, BioPython, BioRuby and other projects' members (2002).

• BioRegistry (Directory)
• Mechanism to specify how and where to retrieve sequence data for each database.
• BioFlat
• Flatfile indexing by using binary tree or BDB(Berkeley DB).
• BioFetch
• Server-client model for getting entry from database via http.
• BioSQL
 Lucas Nussbaum committed Sep 25, 2011 809 
• Schemas to store sequence data to relational databases such as  Lucas Nussbaum committed Jun 15, 2011 810 811 812  MySQL and PostgreSQL, and methods to retrieve entries from the database.
 Lucas Nussbaum committed Sep 25, 2011 813 814 

This tutorial only gives a quick overview of OBDA. Check out the OBDA site for more extensive details.

 Lucas Nussbaum committed Jun 15, 2011 815 816 817 818 819 820 821 822 823 824 825 826 827 828 

BioRegistry

BioRegistry allows for locating retrieval methods and database locations through configuration files. The priorities are

• The file specified with method's parameter
• ~/.bioinformatics/seqdatabase.ini
• /etc/bioinformatics/seqdatabase.ini
• http://www.open-bio.org/registry/seqdatabase.ini

Note that the last locaation refers to www.open-bio.org and is only used when all local configulation files are not available.

In the current BioRuby implementation all local configulation files are read. For databases with the same name settings encountered first are used. This means that if you don't like some settings of a  Lucas Nussbaum committed Sep 25, 2011 829 830 database in the system's global configuration file (/etc/bioinformatics/seqdatabase.ini), you can easily override them by  Lucas Nussbaum committed Jun 15, 2011 831 832 833 834 writing settings to ~/.bioinformatics/seqdatabase.ini.

The syntax of the configuration file is called a stanza format. For example

[DatabaseName] protocol=ProtocolName 
Lucas Nussbaum committed Sep 25, 2011  835     836                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          location=ServerName

You can write a description like the above entry for every database.

 Lucas Nussbaum committed Jun 15, 2011 837 838 839 840 841 842 843 

The database name is a local label for yourself, so you can name it freely and it can differ from the name of the actual databases. In the actual specification of BioRegistry where there are two or more settings for a database of the same name, it is proposed that connection to the database is tried sequentially with the order written in configuration files. However, this has not (yet) been implemented in BioRuby.

 Lucas Nussbaum committed Sep 25, 2011 844 845 

In addition, for some protocols, you must set additional options other than locations (e.g. user name for MySQL). In the BioRegistory  Lucas Nussbaum committed Jun 15, 2011 846 847 848 849 850 851 852 853 854 855 856 857 specification, current available protocols are:

• index-flat
• index-berkeleydb
• biofetch
• biosql
• bsane-corba
• xembl

In BioRuby, you can use index-flat, index-berkleydb, biofetch and biosql. Note that the BioRegistry specification sometimes gets updated and BioRuby does not always follow quickly.

 Lucas Nussbaum committed Sep 25, 2011 858 

Here is an example. It creates a Bio::Registry object and reads the configuration files:

 Lucas Nussbaum committed Jun 15, 2011 859 860 861 862 863 864 865 
reg = Bio::Registry.new  # connects to the database "genbank" serv = reg.get_database('genbank')  # gets entry of the ID entry = serv.get_by_id('AA2CG')
 Lucas Nussbaum committed Sep 25, 2011 866 867 

The variable "serv" is a server object corresponding to the settings written in the configuration files. The class of the object is one of  Lucas Nussbaum committed Jun 15, 2011 868 869 Bio::SQL, Bio::Fetch, and so on. Note that Bio::Registry#get_database("name") returns nil if no database is found.

 Lucas Nussbaum committed Sep 25, 2011 870 871 

After that, you can use the get_by_id method and some specific methods. Please refer to the sections below for more information.

 Lucas Nussbaum committed Jun 15, 2011 872 873 874 

BioFlat

BioFlat is a mechanism to create index files of flat files and to retrieve these entries fast. There are two index types. index-flat is a simple index  Lucas Nussbaum committed Sep 25, 2011 875 performing binary search without using any external libraries of Ruby. index-berkeleydb  Lucas Nussbaum committed Jun 15, 2011 876 uses Berkeley DB for indexing - but requires installing bdb on your computer,  Lucas Nussbaum committed Sep 25, 2011 877 as well as the BDB Ruby package. To create the index itself, you can use br_bioflat.rb command bundled with BioRuby.

 Lucas Nussbaum committed Jun 15, 2011 878 879 
% br_bioflat.rb --makeindex database_name [--format data_format] filename...

The format can be omitted because BioRuby has autodetection. If that  Lucas Nussbaum committed Sep 25, 2011 880 doesn't work, you can try specifying the data format as the name of a BioRuby database class.

 Lucas Nussbaum committed Jun 15, 2011 881 882 

Search and retrieve data from database:

% br_bioflat.rb database_name identifier
 Lucas Nussbaum committed Sep 25, 2011 883 

For example, to create an index of GenBank files gbbct*.seq and get the entry from the database:

 Lucas Nussbaum committed Jun 15, 2011 884 885 886 
% br_bioflat.rb --makeindex my_bctdb --format GenBank gbbct*.seq % br_bioflat.rb my_bctdb A16STM262

If you have Berkeley DB on your system and installed the bdb extension  Lucas Nussbaum committed Sep 25, 2011 887 module of Ruby (see the BDB project page ), you can  Lucas Nussbaum committed Jun 15, 2011 888 889 890 891 892 893 create and search indexes with Berkeley DB - a very fast alternative that uses little computer memory. When creating the index, use the "--makeindex-bdb" option instead of "--makeindex".

% br_bioflat.rb --makeindex-bdb database_name [--format data_format] filename...

BioFetch

Note: this section is an advanced topic
 Lucas Nussbaum committed Sep 25, 2011 894 895 

BioFetch is a database retrieval mechanism via CGI. CGI Parameters, options and error codes are standardized. Client access via  Lucas Nussbaum committed Jun 15, 2011 896 897 http is possible giving the database name, identifiers and format to retrieve entries.

 Lucas Nussbaum committed Sep 25, 2011 898 

The BioRuby project has a BioFetch server at bioruby.org. It uses  Lucas Nussbaum committed Jun 15, 2011 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 GenomeNet's DBGET system as a backend. The source code of the server is in sample/ directory. Currently, there are only two BioFetch servers in the world: bioruby.org and EBI.

Here are some methods to retrieve entries from our BioFetch server.

1. Using a web browser

http://bioruby.org/cgi-bin/biofetch.rb
2. Using the br_biofetch.rb command

% br_biofetch.rb db_name entry_id
3. Directly using Bio::Fetch in a script

serv = Bio::Fetch.new(server_url) entry = serv.fetch(db_name, entry_id)
4. Indirectly using Bio::Fetch via BioRegistry in script

reg = Bio::Registry.new serv = reg.get_database('genbank') entry = serv.get_by_id('AA2CG')
 Lucas Nussbaum committed Sep 25, 2011 916 917 

If you want to use (4), you have to include some settings in seqdatabase.ini. For example:

 Lucas Nussbaum committed Jun 15, 2011 918 919 920 921 922 
[genbank] protocol=biofetch location=http://bioruby.org/cgi-bin/biofetch.rb biodbname=genbank

The combination of BioFetch, Bio::KEGG::GENES and Bio::AAindex1

 Lucas Nussbaum committed Sep 25, 2011 923 924 925 

Bioinformatics is often about gluing things together. Here is an example that gets the bacteriorhodopsin gene (VNG1467G) of the archaea Halobacterium from KEGG GENES database and gets alpha-helix index  Lucas Nussbaum committed Jun 15, 2011 926 data (BURA740101) from the AAindex (Amino acid indices and similarity  Lucas Nussbaum committed Sep 25, 2011 927 matrices) database, and shows the helix score for each 15-aa length  Lucas Nussbaum committed Jun 15, 2011 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 overlapping window.

#!/usr/bin/env ruby  require 'bio'  entry = Bio::Fetch.query('hal', 'VNG1467G') aaseq = Bio::KEGG::GENES.new(entry).aaseq  entry = Bio::Fetch.query('aax1', 'BURA740101') helix = Bio::AAindex1.new(entry).index  position = 1 win_size = 15  aaseq.window_search(win_size) do |subseq|   score = subseq.total(helix)   puts [ position, score ].join("\t")   position += 1 end
 Lucas Nussbaum committed Sep 25, 2011 947 948 

The special method Bio::Fetch.query uses the preset BioFetch server at bioruby.org. (The server internally gets data from GenomeNet.  Lucas Nussbaum committed Jun 15, 2011 949 Because the KEGG/GENES database and AAindex database are not available  Lucas Nussbaum committed Sep 25, 2011 950 from other BioFetch servers, we used the bioruby.org server with  Lucas Nussbaum committed Jun 15, 2011 951 952 Bio::Fetch.query method.)

BioSQL

 Lucas Nussbaum committed Sep 25, 2011 953 954 

BioSQL is a well known schema to store and retrive biological sequences using a RDBMS like PostgreSQL or MySQL: note that SQLite is not supported. First of all, you must install a database engine or have access to a remote one. Then create the schema and populate with the taxonomy. You can follow the Official Guide to accomplish these steps.  Lucas Nussbaum committed Jun 15, 2011 955 956 957 958 959 960 961 Next step is to install these gems:

• ActiveRecord
• CompositePrimaryKeys (Rails doesn't handle by default composite primary keys)
• The layer to comunicate with you preferred RDBMS (postgresql, mysql, jdbcmysql in case you are running JRuby )

You can find ActiveRecord's models in /bioruby/lib/bio/io/biosql

 Lucas Nussbaum committed Sep 25, 2011 962 

When you have your database up and running, you can connect to it like this:

 Lucas Nussbaum committed Jun 15, 2011 963 964 965 966 967 
#!/usr/bin/env ruby  require 'bio'  connection = Bio::SQL.establish_connection({'development'=>{'hostname'=>"YourHostname", 
Lucas Nussbaum committed Sep 25, 2011  968     969     970     971     972     973     974                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  'database'=>"CoolBioSeqDB", 'adapter'=>"jdbcmysql", 'username'=>"YourUser", 'password'=>"YouPassword"       }   }, 'development') 
Lucas Nussbaum committed Jun 15, 2011  975                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
Lucas Nussbaum committed Sep 25, 2011  976     977                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          #The first parameter is the hash contaning the description of the configuration; similar to database.yml in Rails applications, you can declare different environment.  #The second parameter is the environment to use: 'development', 'test', or 'production'. 
Lucas Nussbaum committed Jun 15, 2011  978     979     980     981     982     983     984     985     986     987     988     989     990     991     992     993     994     995                                                                                                                                                                                                                                                                                                                                                                                           #To store a sequence into the database you simply need a biosequence object. biosql_database = Bio::SQL::Biodatabase.find(:first) ff = Bio::GenBank.open("gbvrl1.seq")  ff.each_entry do |gb|   Bio::SQL::Sequence.new(:biosequence=>gb.to_biosequence, :biodatabase=>biosql_database end  #You can list all the entries into every database  Bio::SQL.list_entries  #list databases: Bio::SQL.list_databases  #retriving a generic accession bioseq = Bio::SQL.fetch_accession("YouAccession")  
Lucas Nussbaum committed Sep 25, 2011  996     997                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          #If you use biosequence objects, you will find all its method mapped to BioSQL sequences.  #But you can also access to the models directly: 
Lucas Nussbaum committed Jun 15, 2011  998                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
Lucas Nussbaum committed Sep 25, 2011  999                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  #get the raw sequence associated with your accession 
Lucas Nussbaum committed Jun 15, 2011  1000     1001                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        bioseq.entry.biosequence   
Lucas Nussbaum committed Sep 25, 2011  1002                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 #get the length of your sequence; this is the explicit form of bioseq.length 
Lucas Nussbaum committed Jun 15, 2011  1003     1004                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        bioseq.entry.biosequence.length  
Lucas Nussbaum committed Sep 25, 2011  1005                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 #convert the sequence into GenBank format 
Lucas Nussbaum committed Jun 15, 2011  1006                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 bioseq.to_biosequence.output(:genbank)
 Lucas Nussbaum committed Sep 25, 2011 1007 

BioSQL's schema is not very intuitive for beginners, so spend some time on understanding it. In the end if you know a little bit of Ruby on Rails, everything will go smoothly. You can find information on Annotation here.  Lucas Nussbaum committed Jun 15, 2011 1008 1009 1010 ToDo: add exemaples from George. I remember he did some cool post on BioSQL and Rails.

PhyloXML

PhyloXML is an XML language for saving, analyzing and exchanging data of  Lucas Nussbaum committed Sep 25, 2011 1011 1012 1013 annotated phylogenetic trees. PhyloXML's parser in BioRuby is implemented in Bio::PhyloXML::Parser, and its writer in Bio::PhyloXML::Writer. More information can be found at www.phyloxml.org.

 Lucas Nussbaum committed Jun 15, 2011 1014 

Requirements

 Lucas Nussbaum committed Sep 25, 2011 1015 

In addition to BioRuby, you need the libxml Ruby bindings. To install, execute:

 Lucas Nussbaum committed Jun 15, 2011 1016 
% gem install -r libxml-ruby
 Lucas Nussbaum committed Sep 25, 2011 1017 

 Lucas Nussbaum committed Jun 15, 2011 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 

Parsing a file

require 'bio'  # Create new phyloxml parser phyloxml = Bio::PhyloXML::Parser.open('example.xml')  # Print the names of all trees in the file phyloxml.each do |tree|   puts tree.name end
 Lucas Nussbaum committed Sep 25, 2011 1028 

If there are several trees in the file, you can access the one you wish by specifying its index:

 Lucas Nussbaum committed Jun 15, 2011 1029 
tree = phyloxml[3]
 Lucas Nussbaum committed Sep 25, 2011 1030 

You can use all Bio::Tree methods on the tree, since PhyloXML::Tree inherits from Bio::Tree. For example,

 Lucas Nussbaum committed Jun 15, 2011 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 
tree.leaves.each do |node|   puts node.name end

PhyloXML files can hold additional information besides phylogenies at the end of the file. This info can be accessed through the 'other' array of the parser object.

phyloxml = Bio::PhyloXML::Parser.open('example.xml') while tree = phyloxml.next_tree   # do stuff with trees end   puts phyloxml.other

Writing a file

# Create new phyloxml writer writer = Bio::PhyloXML::Writer.new('tree.xml')  # Write tree to the file tree.xml writer.write(tree1)   # Add another tree to the file writer.write(tree2)

Retrieving data

 Lucas Nussbaum committed Sep 25, 2011 1051 

Here is an example of how to retrieve the scientific name of the clades included in each tree.

 Lucas Nussbaum committed Jun 15, 2011 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 
require 'bio'  phyloxml = Bio::PhyloXML::Parser.open('ncbi_taxonomy_mollusca.xml') phyloxml.each do |tree|   tree.each_node do |node|     print "Scientific name: ", node.taxonomies[0].scientific_name, "\n"   end end

Retrieving 'other' data

require 'bio'  phyloxml = Bio::PhyloXML::Parser.open('phyloxml_examples.xml') while tree = phyloxml.next_tree  #do something with the trees end  p phyloxml.other puts "\n" #=> output is an object representation  #Print in a readable way puts phyloxml.other[0].to_xml, "\n" #=>: # #<align:alignment xmlns:align="http://example.org/align"> #  <seq name="A">acgtcgcggcccgtggaagtcctctcct</seq> #  <seq name="B">aggtcgcggcctgtggaagtcctctcct</seq> #  <seq name="C">taaatcgc--cccgtgg-agtccc-cct</seq> #</align:alignment>  #Once we know whats there, lets output just sequences phyloxml.other[0].children.each do |node|  puts node.value end #=> # #acgtcgcggcccgtggaagtcctctcct #aggtcgcggcctgtggaagtcctctcct #taaatcgc--cccgtgg-agtccc-cct

The BioRuby example programs

 Lucas Nussbaum committed Sep 25, 2011 1092 

Some sample programs are stored in ./samples/ directory. For example, the n2aa.rb program (transforms a nucleic acid sequence into an amino acid sequence) can be run using:

 Lucas Nussbaum committed Jun 15, 2011 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 
./sample/na2aa.rb test/data/fasta/example1.txt

Unit testing and doctests

BioRuby comes with an extensive testing framework with over 1300 tests and 2700 assertions. To run the unit tests:

cd test ruby runner.rb

We have also started with doctest for Ruby. We are porting the examples in this tutorial to doctest - more info upcoming.

See the BioRuby in anger Wiki. A lot of BioRuby's documentation exists in the source code and unit tests. To really dive in you will need the latest source  Lucas Nussbaum committed Sep 25, 2011 1104 code tree. The embedded rdoc documentation for the BioRuby source code can be viewed online at  Lucas Nussbaum committed Jun 15, 2011 1105 1106 <URL:http://bioruby.org/rdoc/>.

BioRuby Shell

 Lucas Nussbaum committed Sep 25, 2011 1107 

The BioRuby shell implementation is located in ./lib/bio/shell. It is very interesting  Lucas Nussbaum committed Jun 15, 2011 1108 as it uses IRB (the Ruby intepreter) which is a powerful environment described in  Lucas Nussbaum committed Sep 25, 2011 1109 Programming Ruby's IRB chapter. IRB commands can be typed directly into the shell, e.g.

 Lucas Nussbaum committed Jun 15, 2011 1110 1111 
bioruby!> IRB.conf[:PROMPT_MODE] ==!> :PROMPT_C
 Lucas Nussbaum committed Sep 25, 2011 1112 

Additionally, you also may want to install the optional Ruby readline support -  Lucas Nussbaum committed Jun 15, 2011 1113 with Debian libreadline-ruby. To edit a previous line you may have to press  Lucas Nussbaum committed Sep 25, 2011 1114 line down (down arrow) first.

 Lucas Nussbaum committed Jun 15, 2011 1115 1116 1117 1118 1119 

Apart from rdoc you may also want to use rtags - which allows jumping around source code by clicking on class and method names.

cd bioruby/lib rtags -R --vi
 Lucas Nussbaum committed Sep 25, 2011 1120 

For a tutorial see here

 Lucas Nussbaum committed Jun 15, 2011 1121 1122 1123 1124 1125 1126 1127 

APPENDIX

KEGG API

Please refer to KEGG_API.rd.ja (English version: <URL:http://www.genome.jp/kegg/soap/doc/keggapi_manual.html> ) and

Ruby Ensembl API

 Lucas Nussbaum committed Sep 25, 2011 1128 

The Ruby Ensembl API is a Ruby API to the Ensembl database. It is NOT currently  Lucas Nussbaum committed Jun 15, 2011 1129 included in the BioRuby archives. To install it, see  Lucas Nussbaum committed Sep 25, 2011 1130 the Ruby-Ensembl Github  Lucas Nussbaum committed Jun 15, 2011 1131 1132 1133 1134 1135 1136 1137 1138 1139 for more information.

Gene Ontology (GO) through the Ruby Ensembl API

Gene Ontologies can be fetched through the Ruby Ensembl API package:

require 'ensembl' Ensembl::Core::DBConnection.connect('drosophila_melanogaster') infile = IO.readlines(ARGV.shift) # reading your comma-separated accession mapping file (one line per mapping) infile.each do |line|   accs = line.split(",")          # Split the comma-sep.entries into an array   drosphila_acc = accs.shift      # the first entry is the Drosophila acc 
Lucas Nussbaum committed Sep 25, 2011  1140                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   mosq_acc = accs.shift           # the second entry is your Mosq. acc 
Lucas Nussbaum committed Jun 15, 2011  1141     1142     1143     1144     1145     1146     1147     1148     1149     1150                                                                                                                                                                                                                                                                                                                                                                                                                                                  gene = Ensembl::Core::Gene.find_by_stable_id(drosophila_acc)   print "#{mosq_acc}"   gene.go_terms.each do |go|      print ",#{go}"   end end

Prints each mosq. accession/uniq identifier and the GO terms from the Drosphila homologues.

Using BioPerl or BioPython from Ruby

At the moment there is no easy way of accessing BioPerl from Ruby. The best way, perhaps, is to create a Perl server that gets accessed through XML/RPC or SOAP.

 Lucas Nussbaum committed Sep 25, 2011 1151 

Installing required external libraries

 Lucas Nussbaum committed Jun 15, 2011 1152 

At this point for using BioRuby no additional libraries are needed, except if  Lucas Nussbaum committed Sep 25, 2011 1153 you are using the Bio::PhyloXML module; then you have to install libxml-ruby.

 Lucas Nussbaum committed Jun 15, 2011 1154 1155 1156 1157 1158 1159 1160 

This may change, so keep an eye on the Bioruby website. Also when a package is missing BioRuby should show an informative message.

At this point installing third party Ruby packages can be a bit painful, as the gem standard for packages evolved late and some still force you to copy things by hand. Therefore read the README's carefully that come with each package.

Installing libxml-ruby

 Lucas Nussbaum committed Sep 25, 2011 1161 

The simplest way is to use the RubyGems packaging system:

 Lucas Nussbaum committed Jun 15, 2011 1162 1163 1164 
gem install -r libxml-ruby

If you get require': no such file to load - mkmf (LoadError) error then do

sudo apt-get install ruby-dev
 Lucas Nussbaum committed Sep 25, 2011 1165 

If you have other problems with installation, then see <URL:http://libxml.rubyforge.org/install.xml>.

 Lucas Nussbaum committed Jun 15, 2011 1166 1167 1168 1169 

Trouble shooting

• Error: in require': no such file to load -- bio (LoadError)
 Lucas Nussbaum committed Sep 25, 2011 1170 

Ruby is failing to find the BioRuby libraries - add it to the RUBYLIB path, or pass  Lucas Nussbaum committed Jun 15, 2011 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 it to the interpeter. For example:

ruby -I\$BIORUBYPATH/lib yourprogram.rb