package Align::NW;

use 5;
use strict;
use vars qw($VERSION);
use base qw(Exporter DynaLoader);

require Exporter;
require DynaLoader;

$VERSION = "1.01";

bootstrap Align::NW $VERSION;


sub new
{
    my($package, $a, $b, $p) = @_;

    my @a = ('', split //, $a);
    my @b = ('', split //, $b);

    my $rows = @a;
    my $cols = @b;

    $a = join '', @a;
    $b = join '', @b;

    my $payoff = new Align::NW::PayoffPtr
	map { $p->{$_} } qw(match mismatch gap_open gap_extend);

    my $matrix = new Align::NW::MatrixPtr $a, $b, $payoff;

    my $nw = { a      => \@a,
	       b      => \@b,
	       rows   => $rows,
	       cols   => $cols,
	       matrix => $matrix,
	       payoff => $payoff };

    bless $nw, $package
}


sub score
{
    my $nw = shift;
    $nw->{matrix}->score;
}

sub align
{
    my $nw = shift;

    $nw->{align} = { a => [],
		     s => [],
		     b => [] };

    my $tail = $nw->max_cell;
	       $nw->align_tail($tail);
    my $head = $nw->align_body($tail);
    	       $nw->align_head($head);
	       $nw->join_align;

    $nw->{score} = $tail->score;
}


sub max(&@)
{
    my $less = shift;
    my $max  = shift;

    &$less($max, $_) and $max = $_  for (@_);

    $max
}


sub max_cell
{
    my $nw     = shift;
    my $matrix = $nw->{matrix};

    my $rows = $nw->{rows};
    my $cols = $nw->{cols};

    my @right  = map { $matrix->cell($_     , $cols-1) } 1..$rows-1;
    my @bottom = map { $matrix->cell($rows-1, $_     ) } 1..$cols-1;

    max { $_[0]->score < $_[1]->score } @right, @bottom
}


sub align_tail
{
    my($nw, $cell) = @_;

    my $a  = $nw->{a};
    my $b  = $nw->{b};

    my $row = $cell->row;
    my $col = $cell->col;

    my @a = @$a[$row+1..$#$a];
    my @s = ();
    my @b = @$b[$col+1..$#$b];
    
    $nw->unshift_align(\@a, \@s, \@b);
}


sub align_body
{
    my($nw, $cell) = @_;

    my $lp = $nw->{lp};
    my $a  = $nw->{a};
    my $b  = $nw->{b};
 
    my(@a, @s, @b);

    for (;;)
    {	
	my $row = $cell->row;
	my $col = $cell->col;
	$row and $col or last;

	my $prev = $cell->prev;

	if ($prev->row < $row and $prev->col < $col)
	{
	    my $a1 = $a->[$row];
	    my $b1 = $b->[$col];
	    unshift @a, $a1;
	    unshift @s, $a1 eq $b1 ? '|' : ' ';
	    unshift @b, $b1;
	}
	elsif ($prev->row < $row)
	{
	    my $gap = $row - $prev->row;
	    unshift @a, @$a[$row-$gap+1..$row];
	    unshift @s, ' ' x $gap;
	    unshift @b, '.' x $gap;
	}
	else
	{
	    my $gap = $col - $prev->col;
	    unshift @a, '.' x $gap;
	    unshift @s, ' ' x $gap;
	    unshift @b, @$b[$col-$gap+1..$col];
	}

	$cell = $prev;
    }

    $nw->unshift_align(\@a, \@s, \@b);
    $cell
}


sub align_head
{
    my($nw, $cell) = @_;

    my $a  = $nw->{a};
    my $b  = $nw->{b};

    my $row = $cell->row;
    my $col = $cell->col;
    my $max = max { $_[0] < $_[1] } $row, $col;

    my @a = (' ' x $col, @$a[1..$row]);
    my @s = (' ' x $max		     );
    my @b = (' ' x $row, @$b[1..$col]);

    $nw->unshift_align(\@a, \@s, \@b);
}


sub unshift_align
{
    my($nw, $a, $s, $b) = @_;
    my $align = $nw->{align};

    unshift @{$align->{a}}, @$a;
    unshift @{$align->{s}}, @$s;
    unshift @{$align->{b}}, @$b;
}


sub join_align
{
    my $nw = shift;
    my $align = $nw->{align};

    for my $key (keys %$align)
    {
	my $x = $nw->{align}{$key};
	$nw->{align}{$key} = join('', @$x);
    }
}    

sub get_score { shift->{score} }
sub get_align { shift->{align} } 


sub print_align
{
    my $nw = shift;
    my $align = $nw->{align};
    my $a = $align->{a};
    my $s = $align->{s};
    my $b = $align->{b};
    my $lineLen = 60;

    $a =~ tr[ -~][^]c;
    $b =~ tr[ -~][^]c;

    for (my $i=0; $i<length($a); $i+=$lineLen)
    {
	print substr($a, $i, $lineLen), "\n";
	print substr($s, $i, $lineLen), "\n";
	print substr($b, $i, $lineLen), "\n";
	print "\n";
    }
}


sub dump_align
{
    my $nw = shift;
    my $align = $nw->{align};

    my $a = $align->{a};
    my $s = $align->{s};
    my $b = $align->{b};

    "$a\n$s\n$b\n"
}


1

__END__

=head1 NAME

Align::NW - Needleman-Wunsch algorithm for optimal global sequence alignment

=head1 SYNOPSIS

    use Align::NW;
  
    $payoff = { match      => $match,
		mismatch   => $mismatch,
		gap_open   => $gap_open,
		gap_extend => $gap_extend };
  
    $nw = new Align::NW $a, $b, $payoff, %options
    $nw->score;
    $nw->align;
  
    $score = $nw->get_score;
    $align = $nw->get_align;

    $nw->print_align;

=head1 DESCRIPTION

C<Align::NW> finds the optimal global alignment of the sequences
C<$a> and C<$b>, subject to the C<$payoff> matrix.

=head2 Algorithm

C<Align::NW> uses the Needleman-Wunsch dynamic programming algorithm.
This algorithm runs in O(a*b*(a+b)), where a and b are the 
lengths of the two sequences to be aligned. 

=head2 Alignments

An alignment of two sequences is represented by three lines.
The top line shows the first sequence,
and the bottom line shows the second sequence.

The middle line has a row of symbols.
The symbol is a vertical bar where ever characters in the two sequences match,
and a space where ever they do not.

Dots may be inserted in either sequence to represent gaps.

For example, the two sequences

    abcdefghajklm
    abbdhijk

could be aligned like this

    abcdefghajklm
    || |   | || 
    abbd...hijk

As shown, there are 6 matches, 2 mismatches, and one gap of length 3.

C<Align::NW> retuns an alignment as a hash

    $align = { a => $a,
	       s => $s,
	       b => $b };

I<$a> and I<$b> are the two sequences. 
I<$s> is the line of symbols.

=head2 The Payoff Matrix

The alignment is scored according to a payoff matrix

    $payoff = { match    => $match,
		mismatch => $mismatch,
		gap_open   	 => $gap_open,
		gap_extend 	 => $gap_extend };

The entries in the matrix are the number of points added to the score 

=over

=item *

for each match

=item *

for each mismatch

=item *

when a gap is opened in either sequence

=item *

for each position that a gap is extended (including the first)

=back

For correct operation, match must be positive, 
and the other entries must be negative.

=head2 Example

Given the payoff matrix

   $payoff = { match      =>  4,
	       mismatch   => -3,
	       gap_open   => -2,
	       gap_extend => -1 };

The sequences

    abcdefghajklm
    abbdhijk

are aligned and scored like this

                a b c d e f g h a j k l m
                | |   |       |   | | 
                a b b d . . . h i j k

    match       4 4   4       4   4 4  
    mismatch       -3          -3
    gap_open           -2
    gap_extend         -1-1-1

for a total score of 24-6-2-3 = 15. 
The algorithm guarantees that no other alignment of these two sequences
has a higher score under this payoff matrix.

=head1 METHODS

=over 4

=item I<$nw> = C<new> C<Align::NW> I<$a>, I<$b>, I<$payoff>, I<%options>

Creates and returns a new C<Align::NW> object.
I<$a> and I<$b> are the sequences to be aligned.
I<$payoff> is the payoff matrix, described above.
Additional options maybe passed in the I<%options> hash;
see L</OPTIONS> for details.

=item I<$nw>->C<score>

Fills in the score matrix for I<$nw>.
This is the O(a*b*(a+b)) operation.

=item I<$nw>->C<align>

Backtracks through the score matrix and generates an alignment for 
the two sequences.
C<score> must be called before C<align>.

=item I<$score> = I<$nw>->C<get_score>

Returns the score of the alignment.
C<score> must be called before C<get_score>.

=item I<$align> = I<$nw>->C<get_align>

Returns the alignment of the two sequences, as described above in 
L</Alignments>.
C<align> must be called before C<get_align>.

=item I<$nw>->C<print_align>

Pretty prints the alignment to STDOUT.
C<align> must be called before C<print_align>.

=back

=head1 OPTIONS

Options may be passed to C<new> in the C<%options> hash.
The following options are defined.

=over 4

=item B<-v>

Verbose output.
Prints some dots to STDERR.
Useful for monitoring the progress of large alignments.

=back

=head1 SEE ALSO

=over 4

=item *

Needleman, S.B. and Wunsch, C.D. 1970. "A general method applicable to
the search for similarities in the amino acid sequences of two
proteins" I<Journal of Molecular Biology>. 48: 443-453.

=item *

Smith, T.F. and Waterman, M.S. 1981. "Identification of common
molecular subsequences" I<Journal of Molecular Biology>. 147: 195-197

=back

There are usually some some tutorials on Needleman-Wunsch and
Smith-Waterman alignment floating around on the web. I used to provide
links to some, but they kept going 404. If you Google around a bit you
can probably find a current one.

=head1 ACKNOWLEDGMENTS

=over 4

=item *

Andreas Doms <ad11@inf.tu-dresden.de>

=back

=head1 AUTHOR

Steven McDougall, <swmcd@world.std.com>

=head1 COPYRIGHT

Copyright 1999-2003 by Steven McDougall. This module is free
software; you can redistribute it and/or modify it under the same
terms as Perl itself.

=cut