NAME
Bio::Tools::dpAlign - Perl extension to do pairwise dynamic programming sequence alignment
SYNOPSIS
use Bio::Tools::dpAlign;
use Bio::SeqIO;
use Bio::SimpleAlign;
use Bio::AlignIO;
$seq1 = new Bio::SeqIO(-file => $ARGV[0], -format => 'fasta');
$seq2 = new Bio::SeqIO(-file => $ARGV[1], -format => 'fasta');
# create a dpAlign object
$factory = new dpAlign(-match => 3,
-mismatch => -1,
-gap => 3,
-ext => 1,
-alg => Bio::Tools::dpAlign::DPALIGN_LOCAL_MILLER_MYERS);
# actually do the alignment
$out = $factory->pairwise_alignment($seq1->next_seq, $seq2->next_seq);
$alnout = new Bio::AlignIO(-format => 'pfam', -fh => \*STDOUT);
$alnout->write_aln($out);
# To do protein alignment, set the sequence type to protein
# currently all protein alignments are using BLOSUM62 matrix
# the gap opening cost is 10 and gap extension is 2. These
# values are from ssearch. They won't be changed even though
# you set other values for now. Also DPALIGN_LOCAL_GREEN is not
# supported for protein in this version.
$seq1->alphabet('protein');
$seq2->alphabet('protein');
$out = $factory->pairwise_alignment($seq1->next_seq, $seq2->next_seq);
$alnout->write_aln($out);
# use the factory to make some output
$factory->align_and_show($seq1, $seq2, STDOUT);DESCRIPTION
Dynamic Programming approach is considered to be the most
sensitive way to align two biological sequences. There are
currently three major types of dynamic programming algorithms:
Global Alignment, Local Alignment and Ends-free Alignment.
Global Alignment compares two sequences in their entirety.
By inserting gaps in the two sequences, it aligns two
sequences to minimize the edit distance as defined by the
gap cost function and the substitution matrix. Global Alignment
is generally applied to two sequences that are very similar
in length and content.
Local Alignment instead attempts to find out the subsequences
that has the minimal edit distance among all possible subsequences.
It is good for sequences that has a stretch of subsequences
that are similar to each other.
Ends-free Alignment is a special case of Global Alignment. There
are no gap penalty imposed for the gaps that extended from
the end points of two sequences. Therefore it will be a good
application when you think one sequence is contained by the
other or when you think two sequences overlap each other.
Dynamic Programming was first introduced by Needleman-Wunsch (1970)
to globally align two sequences. The idea of local alignment
was later introduced by Smith-Waterman (1981). Gotoh (1982)
improved both algorithms by introducing auxillary arrays that
reduced the time complexity of the algorithms to O(m*n).
Miller-Myers (1988) exploits the divide-and-conquer idea
introduced by Hirschberg (1975) to solve the affine gap cost
dynamic programming using only linear space. At the time of
this writing, it is accepted that Miller-Myers is the fastest
single CPU implementation and using the least memory that is
truly equivalent to original algorithm introduced by
Needleman-Wunsch. According to Aaron Mackey, Phil Green's SWAT
implemention introduced a heuristic that does not consider
paths throught the matrix where the score would be less than
the gap opening penalty, yielding a 1.5-2X speedup on most
comparisons. to skip the calculation of some cells. However,
his approach is only good for calculating the minimum edit
distance and find out the corresponding subsequences (aka
search phase). Bill Pearson's popular dynamic programming
alignment program SSEARCH uses Phil Green's algorithm to
find the subsequences and then Miller-Myers's algorithm to
find the actual alignment. (aka alignment phase)
The current implementation supports local alignment of
either DNA sequences or protein sequences. It allows you
to specify either the Phil Green (DPALIGN_LOCAL_GREEN)
or Miller-Myers (DPALIGN_LOCAL_MILLER_MYERS). For DNA
alignment, you can specify the scores for match, mismatch,
gap opening cost and gap extension cost. For protein
alignment, it is using BLOSUM62 by default. Currently the
substitution matrix is not configurable.DEPENDENCIES
This package comes with the main bioperl distribution. You
also need to install the lastest bioperl-ext package which
contains the XS code that implements the algorithms. This
package won't work if you haven't compiled the bioperl-ext
package.TO-DO
1) Allow custom substitution matrix.
2) Support Global Alignment.
3) Support Ends-free Alignment.
4) Include a score only calculation based on Phil Green's
algorithm. The code will be borrowed from do_work in
ssearch.
5) Support IUPAC code for DNA sequenceFEEDBACK
Mailing Lists
User feedback is an integral part of the evolution of this and other Bioperl modules. Send your comments and suggestions preferably to one of the Bioperl mailing lists. Your participation is much appreciated.
bioperl-l@bioperl.org - General discussion
http://bioperl.org/MailList.shtml - About the mailing listsReporting Bugs
Report bugs to the Bioperl bug tracking system to help us keep track the bugs and their resolution. Bug reports can be submitted via email or the web:
bioperl-bugs@bio.perl.org
http://bugzilla.bioperl.org/AUTHOR
This implementation was written by Yee Man Chan (ymc@yahoo.com).
Copyright (c) 2003 Yee Man Chan. All rights reserved. This program
is free software; you can redistribute it and/or modify it under
the same terms as Perl itself. Special thanks to Aaron Mackey
and WIlliam Pearson for the helpful discussions. [The portion
of code inside pgreen subdirectory was borrowed from ssearch. It
should be distributed in the same terms as ssearch.]pairwise_alignment
Title : pairwise_alignment
Usage : $aln = $factory->pairwise_alignment($seq1,$seq2)
Function: Makes a SimpleAlign object from two sequences
Returns : A SimpleAlign object
Args :align_and_show
Title : align_and_show
Usage : $factory->align_and_show($seq1,$seq2,STDOUT)match
Title : match
Usage : $match = $factory->match() #get
: $factory->match($value) #set
Function : the set get for the match score
Example :
Returns : match value
Arguments : new valuemismatch
Title : mismatch
Usage : $mismatch = $factory->mismatch() #get
: $factory->mismatch($value) #set
Function : the set get for the mismatch penalty
Example :
Returns : mismatch value
Arguments : new valuegap
Title : gap
Usage : $gap = $factory->gap() #get
: $factory->gap($value) #set
Function : the set get for the gap penalty
Example :
Returns : gap value
Arguments : new valueext
Title : ext
Usage : $ext = $factory->ext() #get
: $factory->ext($value) #set
Function : the set get for the ext penalty
Example :
Returns : ext value
Arguments : new valuealg
Title : alg
Usage : $alg = $factory->alg() #get
: $factory->alg($value) #set
Function : the set get for the algorithm
Example :
Returns : alg value
Arguments : new value