1 #summary Tutorial for multiple sequence alignments and phylogenetic methods in BioRuby -- under development!
9 Tutorial for multiple sequence alignments and phylogenetic methods in [http://bioruby.open-bio.org/ BioRuby].
11 Eventually, this is expected to be placed on the official !BioRuby page.
13 Author: [https://sites.google.com/site/cmzmasek/ Christian Zmasek], Sanford-Burnham Medical Research Institute
16 Copyright (C) 2011 Christian M Zmasek. All rights reserved.
19 = Multiple Sequence Alignment =
22 == Multiple Sequence Alignment Input and Output ==
24 === Reading in a Multiple Sequence Alignment from a File ===
26 This automatically determines the format
32 ff = Bio::FlatFile.auto('bcl2.fasta')
33 ff.each_entry do |entry|
35 puts entry.entry_id # prints the identifier of the entry
36 puts entry.definition # prints the definition of the entry
37 puts entry.seq # prints the sequence data of the entry
40 # Creates a multiple sequence alignment (possibly unaligned) named
41 # 'seqs' from array 'seq_ary'.
42 seqs = Bio::Alignment.new(seq_ary)
43 seqs.each { |seq| puts seq.to_s }
45 # Writes multiple sequence alignment (possibly unaligned) 'seqs'
46 # to a file in PHYLIP format.
47 File.open('out0.phylip', 'w') do |f|
48 f.write(seqs.output(:phylip))
51 # Writes multiple sequence alignment (possibly unaligned) 'seqs'
52 # to a file in FASTA format.
53 File.open('out0.fasta', 'w') do |f|
54 f.write(seqs.output(:fasta))
59 ==== ClustalW Format ====
61 The following example shows how to read in a *ClustalW*-formatted multiple sequence alignment.
67 # Reads in a ClustalW-formatted multiple sequence alignment
68 # from a file named "infile_clustalw.aln" and stores it in 'report'.
69 report = Bio::ClustalW::Report.new(File.read('infile_clustalw.aln'))
71 # Accesses the actual alignment.
72 msa = report.alignment
74 # Goes through all sequences in 'msa' and prints the
75 # actual molecular sequence.
81 ==== FASTA Format ====
83 The following example shows how to read in a *FASTA*-formatted multiple sequence file. (_This seems a little clumsy, I wonder if there is a more direct way, avoiding the creation of an array.)
88 # Reads in a FASTA-formatted multiple sequence alignment (which does
89 # not have to be aligned, though) and stores its sequences in
92 fasta_seqs = Bio::Alignment::MultiFastaFormat.new(File.open('infile.fasta').read)
93 fasta_seqs.entries.each do |seq|
97 # Creates a multiple sequence alignment (possibly unaligned) named
98 # 'seqs' from array 'seq_ary'.
99 seqs = Bio::Alignment.new(seq_ary)
101 # Prints each sequence to the console.
102 seqs.each { |seq| puts seq.to_s }
104 # Writes multiple sequence alignment (possibly unaligned) 'seqs'
105 # to a file in PHYLIP format.
106 File.open('outfile.phylip', 'w') do |f|
107 f.write(seqs.output(:phylip))
111 Relevant API documentation:
113 * [http://bioruby.open-bio.org/rdoc/classes/Bio/ClustalW/Report.html Bio::ClustalW::Report]
114 * [http://bioruby.open-bio.org/rdoc/classes/Bio/Alignment.html Bio::Alignment]
115 * [http://bioruby.open-bio.org/rdoc/classes/Bio/Sequence.html Bio::Sequence]
117 === Creating a Multiple Sequence Alignment ===
120 === Creating a Multiple Sequence Alignment from a Database ===
124 === Writing a Multiple Sequence Alignment to a File ===
127 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').
133 # Creates a new file named "outfile.fasta" and writes
134 # multiple sequence alignment 'msa' to it in fasta format.
135 File.open('outfile.fasta', 'w') do |f|
136 f.write(msa.output(:fasta))
140 ==== Setting the Output Format ====
142 The following symbols determine the output format:
144 * `:clustal` for ClustalW
146 * `:phylip` for PHYLIP interleaved (will truncate sequence names to no more than 10 characters)
147 * `:phylipnon` for PHYLIP non-interleaved (will truncate sequence names to no more than 10 characters)
149 * `:molphy` for Molphy
152 For example, the following writes in PHYLIP's non-interleaved format:
155 f.write(align.output(:phylipnon))
159 === Formatting of Individual Sequences ===
161 !BioRuby can format molecular sequences in a variety of formats.
162 Individual sequences can be formatted to (e.g.) Genbank format as shown in the following examples.
164 For Sequence objects:
166 seq.to_seq.output(:genbank)
169 For Bio::!FlatFile entries:
171 entry.to_biosequence.output(:genbank)
174 The following symbols determine the output format:
175 * `:genbank` for Genbank
178 * `:fasta_ncbi` for NCBI-type FASTA
179 * `:raw` for raw sequence
180 * `:fastq` for FASTQ (includes quality scores)
181 * `:fastq_sanger` for Sanger-type FASTQ
182 * `:fastq_solexa` for Solexa-type FASTQ
183 * `:fastq_illumina` for Illumina-type FASTQ
185 == Calculating Multiple Sequence Alignments ==
187 !BioRuby can be used to execute a variety of multiple sequence alignment
188 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]).
189 In the following, examples for using the MAFFT and Muscle are shown.
194 The following example uses the MAFFT program to align four sequences
195 and then prints the result to the screen.
196 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.
202 # 'seqs' is either an array of sequences or a multiple sequence
203 # alignment. In general this is read in from a file as described in ?.
204 # For the purpose of this tutorial, it is generated in code.
205 seqs = ['MFQIPEFEPSEQEDSSSAER',
206 'MGTPKQPSLAPAHALGLRKS',
211 # Calculates the alignment using the MAFFT program on the local
212 # machine with options '--maxiterate 1000 --localpair'
213 # and stores the result in 'report'.
214 options = ['--maxiterate', '1000', '--localpair']
215 mafft = Bio::MAFFT.new('path/to/mafft', options)
216 report = mafft.query_align(seqs)
218 # Accesses the actual alignment.
219 align = report.alignment
221 # Prints each sequence to the console.
222 align.each { |s| puts s.to_s }
228 * Katoh, Toh (2008) "Recent developments in the MAFFT multiple sequence alignment program" Briefings in Bioinformatics 9:286-298
230 * Katoh, Toh 2010 (2010) "Parallelization of the MAFFT multiple sequence alignment program" Bioinformatics 26:1899-1900
240 # 'seqs' is either an array of sequences or a multiple sequence
241 # alignment. In general this is read in from a file as described in ?.
242 # For the purpose of this tutorial, it is generated in code.
243 seqs = ['MFQIPEFEPSEQEDSSSAER',
244 'MGTPKQPSLAPAHALGLRKS',
248 # Calculates the alignment using the Muscle program on the local
249 # machine with options '-quiet -maxiters 64'
250 # and stores the result in 'report'.
251 options = ['-quiet', '-maxiters', '64']
252 muscle = Bio::Muscle.new('path/to/muscle', options)
253 report = muscle.query_align(seqs)
255 # Accesses the actual alignment.
256 align = report.alignment
258 # Prints each sequence to the console.
259 align.each { |s| puts s.to_s }
265 * Edgar, R.C. (2004) "MUSCLE: multiple sequence alignment with high accuracy and high throughput" Nucleic Acids Res 32(5):1792-1797
267 === Other Programs ===
269 _need more detail here..._
271 [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.
275 == Manipulating Multiple Sequence Alignments ==
277 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.
291 = Phylogenetic Trees =
294 == Phylogenetic Tree Input and Output ==
297 === Reading in of Phylogenetic Trees ===
301 ====Newick or New Hampshire Format====
311 ====phyloXML Format====
313 Partially copied from [https://www.nescent.org/wg_phyloinformatics/BioRuby_PhyloXML_HowTo_documentation Diana Jaunzeikare's documentation].
315 In addition to !BioRuby, a libxml Ruby binding is also required. This can be installed with the following command:
318 % gem install -r libxml-ruby
321 This example reads file "example.xml" and stores its [http://www.phyloxml.org/ phyloXML]-formatted trees in variable 'trees'.
327 # This creates new phyloXML parser.
328 trees = Bio::PhyloXML::Parser.new('example.xml')
330 # This prints the names of all trees in the file.
335 # If there are several trees in the file, you can access the one you wish via index.
351 === Writing of Phylogenetic Trees ===
353 ====Newick or New Hampshire Format====
363 ====phyloXML Format====
365 Partially copied from [https://www.nescent.org/wg_phyloinformatics/BioRuby_PhyloXML_HowTo_documentation Diana Jaunzeikare's documentation].
367 In addition to !BioRuby, a libxml Ruby binding is also required. This can be installed with the following command:
370 % gem install -r libxml-ruby
373 This example writes trees 'tree1' and 'tree2' to file "tree.xml" in [http://www.phyloxml.org/ phyloXML] format.
379 # this creates new phyloXML writer.
380 writer = Bio::PhyloXML::Writer.new('tree.xml')
382 # Writes tree to the file "tree.xml".
385 # Adds another tree to the file.
403 = Phylogenetic Inference =
405 _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..._
407 == Optimality Criteria Based on Character Data ==
409 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.
411 === Maximum Likelihood ===
434 === Maximum Parsimony ===
436 Currently no direct support in !BioRuby.
439 === Bayesian Inference ===
441 Currently no direct support in !BioRuby.
444 == Pairwise Distance Based Methods ==
446 === Pairwise Sequence Distance Estimation ===
457 === Optimality Criteria Based on Pairwise Distances ===
460 ==== Minimal Evolution: FastME ====
470 === Algorithmic Methods Based on Pairwise Distances ===
472 ==== Neighbor Joining and Related Methods ====
488 == Support Calculation? ==
490 === Bootstrap Resampling? ===
495 = Analyzing Phylogenetic Trees =
500 == Gene Duplication Inference ==
502 _need to further test and then import GSoC 'SDI' work..._
510 = Putting It All Together =
512 Example of a small "pipeline"-type program running a mininal phyogenetic analysis: starting with a set of sequences and ending with a phylogenetic tree.