ab6faf9d1c6fc050e9b264ed01fd4959810a9161
[jalview.git] / wiki / PhyloBioRuby.wiki
1 #summary Tutorial for multiple sequence alignments and phylogenetic methods in BioRuby -- under development!
2
3
4
5 = Introduction =
6
7 Under development!
8
9 Tutorial for multiple sequence alignments and phylogenetic methods in [http://bioruby.open-bio.org/ BioRuby].
10
11 Eventually, this is expected to be placed on the official !BioRuby page.
12
13 Author: [http://www.cmzmasek.net/ Christian M Zmasek], Sanford-Burnham Medical Research Institute
14
15  
16 Copyright (C) 2011 Christian M Zmasek. All rights reserved.
17
18
19 = Multiple Sequence Alignment =
20
21
22 == Multiple Sequence Alignment Input and Output ==
23
24 === Reading in a Multiple Sequence Alignment from a File ===
25
26 The following example shows how to read in a *ClustalW*-formatted multiple sequence alignment.
27
28 {{{
29 #!/usr/bin/env ruby
30 require 'bio'
31
32 # Reads in a ClustalW-formatted multiple sequence alignment
33 # from a file named "infile_clustalw.aln" and stores it in 'report'.
34 report = Bio::ClustalW::Report.new(File.read('infile_clustalw.aln'))
35
36 # Accesses the actual alignment.
37 msa = report.alignment
38
39 # Goes through all sequences in 'msa' and prints the
40 # actual molecular sequence.
41 msa.each do |entry|
42   puts entry.seq
43 end
44 }}}
45
46 Blah
47 {{{
48 # Reads in a Fasta-formatted multiple sequence alignment (which does
49 # not have to be aligned, though) and stores its sequences in
50 # array 'seq_ary'.
51 seq_ary = Array.new
52 fasta_seqs = Bio::Alignment::MultiFastaFormat.new(File.open('bcl2.fasta').read)
53 fasta_seqs.entries.each do |seq|
54   seq_ary.push( seq )
55 end
56
57 # Creates a multiple sequence alignment (possibly unaligned) named
58 # 'seqs' from array 'seq_ary'.
59 seqs = Bio::Alignment.new( seq_ary )
60 seqs.each { |seq| puts seq.to_s }
61
62
63 puts seqs.consensus
64
65 # Writes multiple sequence alignment (possibly unaligned) 'seqs'
66 # to a file in phylip format.
67 File.open('out1.phylip', 'w') do |f|
68   f.write(seqs.output(:phylip))
69 end
70 }}}
71
72 Relevant API documentation:
73
74  * [http://bioruby.open-bio.org/rdoc/classes/Bio/ClustalW/Report.html Bio::ClustalW::Report]
75  * [http://bioruby.open-bio.org/rdoc/classes/Bio/Alignment.html Bio::Alignment]
76  * [http://bioruby.open-bio.org/rdoc/classes/Bio/Sequence.html Bio::Sequence]
77
78 === Writing a Multiple Sequence Alignment to a File ===
79
80
81 The following example shows how to write a multiple sequence alignment in *FASTA*-format. It first creates a file named "outfile.fasta" for writing ('w') and then writes the multiple sequence alignment referred to by variable 'msa' to it in FASTA-format (':fasta').
82
83 {{{
84 #!/usr/bin/env ruby
85 require 'bio'
86
87 # Creates a new file named "outfile.fasta" and writes
88 # multiple sequence alignment 'msa' to it in fasta format.
89 File.open('outfile.fasta', 'w') do |f|
90   f.write(msa.output(:fasta))
91 end
92 }}}
93
94 ==== Setting the Output Format ====
95
96 The following symbols determine the output format:
97
98   * `:clustal` for ClustalW
99   * `:fasta` for FASTA
100   * `:phylip` for PHYLIP interleaved (will truncate sequence names to no more than 10 characters)
101   * `:phylipnon` for PHYLIP non-interleaved (will truncate sequence names to no more than 10 characters)
102   * `:msf` for MSF
103   * `:molphy` for Molphy
104
105
106 For example, the following writes in PHYLIP's non-interleaved format:
107
108 {{{
109 f.write(align.output(:phylipnon))
110 }}}
111
112
113 === Formatting of Individual Sequences ===
114
115 !BioRuby can format molecular sequences in a variety of formats.
116 Individual sequences can be formatted to (e.g.) Genbank format as shown in the following examples.
117
118 For Sequence objects:
119 {{{
120 seq.to_seq.output(:genbank)
121 }}}
122
123 For Bio::!FlatFile entries:
124 {{{
125 entry.to_biosequence.output(:genbank)
126 }}}
127
128 The following symbols determine the output format:
129   * `:genbank` for Genbank
130   * `:embl` for EMBL
131   * `:fasta` for FASTA
132   * `:fasta_ncbi` for NCBI-type FASTA
133   * `:raw` for raw sequence
134   * `:fastq` for FASTQ (includes quality scores)
135   * `:fastq_sanger` for Sanger-type FASTQ 
136   * `:fastq_solexa` for Solexa-type FASTQ 
137   * `:fastq_illumina` for Illumina-type FASTQ 
138
139 == Calculating Multiple Sequence Alignments ==
140
141 !BioRuby can be used to execute a variety of multiple sequence alignment
142 programs (such as [http://mafft.cbrc.jp/alignment/software/ MAFFT], [http://probcons.stanford.edu/ Probcons], [http://www.clustal.org/ ClustalW], [http://www.drive5.com/muscle/ Muscle], and [http://www.tcoffee.org/Projects_home_page/t_coffee_home_page.html T-Coffee]). 
143 In the following, examples for using the MAFFT and Muscle are shown.
144
145
146 === MAFFT ===
147
148 The following example uses the MAFFT program to align four sequences
149 and then prints the result to the screen.
150 Please note that if the path to the MAFFT executable is properly set `mafft=Bio::MAFFT.new(options)` can be used instead of explicitly indicating the path as in the example. 
151
152 {{{
153 #!/usr/bin/env ruby
154 require 'bio'
155
156 # 'seqs' is either an array of sequences or a multiple sequence 
157 # alignment. In general this is read in from a file as described in ?.
158 # For the purpose of this tutorial, it is generated in code.
159 seqs = ['MFQIPEFEPSEQEDSSSAER',
160         'MGTPKQPSLAPAHALGLRKS',
161         'PKQPSLAPAHALGLRKS',
162         'MCSTSGCDLE'] 
163
164
165 # Calculates the alignment using the MAFFT program on the local
166 # machine with options '--maxiterate 1000 --localpair'
167 # and stores the result in 'report'.
168 options = ['--maxiterate', '1000', '--localpair']
169 mafft = Bio::MAFFT.new('path/to/mafft', options)
170 report = mafft.query_align(seqs)
171
172 # Accesses the actual alignment.
173 align = report.alignment
174
175 # Prints each sequence to the console.
176 align.each { |s| puts s.to_s }
177
178 }}}
179
180 References:
181
182  * Katoh, Toh (2008) "Recent developments in the MAFFT multiple sequence alignment program" Briefings in Bioinformatics 9:286-298 
183
184  * Katoh, Toh 2010 (2010) "Parallelization of the MAFFT multiple sequence alignment program" Bioinformatics 26:1899-1900 
185
186
187
188 === Muscle ===
189
190 {{{
191 #!/usr/bin/env ruby
192 require 'bio'
193
194 # 'seqs' is either an array of sequences or a multiple sequence 
195 # alignment. In general this is read in from a file as described in ?.
196 # For the purpose of this tutorial, it is generated in code.
197 seqs = ['MFQIPEFEPSEQEDSSSAER',
198         'MGTPKQPSLAPAHALGLRKS',
199         'PKQPSLAPAHALGLRKS',
200         'MCSTSGCDLE']  
201
202 # Calculates the alignment using the Muscle program on the local
203 # machine with options '-quiet -maxiters 64'
204 # and stores the result in 'report'.
205 options = ['-quiet', '-maxiters', '64']
206 muscle = Bio::Muscle.new('path/to/muscle', options)
207 report = muscle.query_align(seqs)
208
209 # Accesses the actual alignment.
210 align = report.alignment
211
212 # Prints each sequence to the console.
213 align.each { |s| puts s.to_s }
214
215 }}}
216
217 References:
218
219  * Edgar, R.C. (2004) "MUSCLE: multiple sequence alignment with high accuracy and high throughput" Nucleic Acids Res 32(5):1792-1797
220
221 === Other Programs ===
222
223 _need more detail here..._
224
225 [http://probcons.stanford.edu/ Probcons], [http://www.clustal.org/ ClustalW], and [http://www.tcoffee.org/Projects_home_page/t_coffee_home_page.html T-Coffee] can be used in the same manner as the programs above. 
226
227
228
229 == Manipulating Multiple Sequence Alignments ==
230
231 Oftentimes, multiple sequence to be used for phylogenetic inference are 'cleaned up' in some manner. For instance, some researchers prefer to delete columns with more than 50% gaps. The following code is an example of how to do that in !BioRuby.
232
233
234 _... to be done_
235
236 {{{
237 #!/usr/bin/env ruby
238 require 'bio'
239
240 }}}
241
242
243 ----
244
245 = Phylogenetic Trees =
246
247
248 == Phylogenetic Tree Input and Output ==
249
250
251 === Reading in of Phylogenetic Trees ===
252
253
254
255 ====Newick or New Hampshire Format====
256
257 _... to be done_
258
259 {{{
260 #!/usr/bin/env ruby
261 require 'bio'
262
263 }}}
264
265 ====phyloXML Format====
266
267 Partially copied from [https://www.nescent.org/wg_phyloinformatics/BioRuby_PhyloXML_HowTo_documentation Diana Jaunzeikare's documentation].
268
269 In addition to !BioRuby, a libxml Ruby binding is also required. This can be installed with the following command: 
270
271 {{{
272 % gem install -r libxml-ruby
273 }}}
274
275 This example reads file "example.xml" and stores its [http://www.phyloxml.org/ phyloXML]-formatted trees in variable 'trees'.
276
277 {{{
278 #!/usr/bin/env ruby
279 require 'bio'
280
281 # This creates new phyloXML parser.
282 trees = Bio::PhyloXML::Parser.new('example.xml')
283
284 # This prints the names of all trees in the file.
285 trees.each do |tree|
286   puts tree.name
287 end
288
289 # If there are several trees in the file, you can access the one you wish via index.
290 tree = trees[3]
291
292 }}}
293
294
295 ====Nexus  Format====
296
297 _... to be done_
298
299 {{{
300 #!/usr/bin/env ruby
301 require 'bio'
302
303 }}}
304
305 === Writing of Phylogenetic Trees ===
306
307 ====Newick or New Hampshire Format====
308
309 _... to be done_
310
311 {{{
312 #!/usr/bin/env ruby
313 require 'bio'
314
315 }}}
316
317 ====phyloXML Format====
318
319 Partially copied from [https://www.nescent.org/wg_phyloinformatics/BioRuby_PhyloXML_HowTo_documentation Diana Jaunzeikare's documentation].
320
321 In addition to !BioRuby, a libxml Ruby binding is also required. This can be installed with the following command: 
322
323 {{{
324 % gem install -r libxml-ruby
325 }}}
326
327 This example writes trees 'tree1' and 'tree2' to file "tree.xml" in [http://www.phyloxml.org/ phyloXML] format.
328
329 {{{
330 #!/usr/bin/env ruby
331 require 'bio'
332
333 # this creates new phyloXML writer.
334 writer = Bio::PhyloXML::Writer.new('tree.xml')
335
336 # Writes tree to the file "tree.xml".
337 writer.write(tree1)
338
339 # Adds another tree to the file.
340 writer.write(tree2)
341
342
343 }}}
344
345
346 ====Nexus  Format====
347
348 _... to be done_
349
350 {{{
351 #!/usr/bin/env ruby
352 require 'bio'
353
354 }}}
355
356
357 = Phylogenetic Inference =
358
359 _Currently !BioRuby does not contain wrappers for phylogenetic inference programs, thus I am progress of writing a RAxML wrapper followed by a wrapper for FastME..._
360
361 == Optimality Criteria Based on Character Data ==
362
363 Character data based methods work directly on molecular sequences and thus do not require the calculation of pairwise distances but tend to be time consuming and sensitive to errors in the multiple sequence alignment.
364
365 === Maximum Likelihood ===
366
367 ==== RAxML ====
368
369 _... to be done_
370
371 {{{
372 #!/usr/bin/env ruby
373 require 'bio'
374
375 }}}
376
377
378 ==== PhyML ====
379
380 _... to be done_
381
382 {{{
383 #!/usr/bin/env ruby
384 require 'bio'
385
386 }}}
387
388 === Maximum Parsimony ===
389
390 Currently no direct support in !BioRuby.
391
392
393 === Bayesian Inference ===
394
395 Currently no direct support in !BioRuby.
396
397
398 == Pairwise Distance Based Methods ==
399
400 === Pairwise Sequence Distance Estimation ===
401
402 _... to be done_
403
404 {{{
405 #!/usr/bin/env ruby
406 require 'bio'
407
408 }}}
409
410
411 === Optimality Criteria Based on Pairwise Distances ===
412
413
414 ==== Minimal Evolution: FastME ====
415
416 _... to be done_
417
418 {{{
419 #!/usr/bin/env ruby
420 require 'bio'
421
422 }}}
423
424 === Algorithmic Methods Based on Pairwise Distances ===
425
426 ==== Neighbor Joining and Related Methods ====
427
428 _... to be done_
429
430 {{{
431 #!/usr/bin/env ruby
432 require 'bio'
433
434 }}}
435
436
437
438
439
440
441
442 == Support Calculation? ==
443
444 === Bootstrap Resampling? ===
445
446
447 ----
448
449 = Analyzing Phylogenetic Trees =
450
451 == PAML ==
452
453
454 == Gene Duplication Inference ==
455
456 _need to further test and then import GSoC 'SDI' work..._
457
458
459 == Others? ==
460
461
462 ----
463
464 = Putting It All Together =
465
466 Example of a small "pipeline"-type program running a mininal phyogenetic analysis: starting with a set of sequences and ending with a phylogenetic tree.
467