view release on metacpan or  search on metacpan

lib/Algorithm/NeedlemanWunsch.pm  view on Meta::CPAN

	while ($j <= $n) {
	    $D->[0]->[$j] = $j * $self->{gap_penalty};
	    ++$j;
	}
    } else {
        $j = 1;
	while ($j <= $n) {
	    $D->[0]->[$j] = 0;
	    ++$j;
        }
    }

    my $i = 1;
    while ($i <= $m) {
        $A->[$i]->[0] = $from_up;
	++$i;
    }

    $i = 1;
    while ($i <= $m) {
	$D->[$i]->[0] = $i * $self->{gap_penalty};
	++$i;
    }

    $i = 1;
    while ($i <= $m) {
	$j = 1;
	while ($j <= $n) {
	    my @scores = map { &$_($i, $j); } @subproblems;
	    my $d = max(@scores);

	    my $a = 0;
	    my $from = 1;
	    my $k = 0;
	    while ($k < scalar(@scores)) {
		if ($scores[$k] == $d) {
		  $a |= $from;
		}

		$from *= 2;
		++$k;
	    }

	    $A->[$i]->[$j] = $a;
	    $D->[$i]->[$j] = $d;

	    # my $x = join ', ', @scores;
	    # warn "$i, $j: $x -> ", $D->[$i]->[$j], "\n";
	    ++$j;
	}

	++$i;
    }

    $i = $m;
    $j = $n;
    while (($i > 0) || ($j > 0)) {
        my $a = $A->[$i]->[$j];
	my @alt;
	if ($a & $from_diag) {
	    die "internal error" unless ($i > 0) && ($j > 0);
	    push @alt, [ $i - 1, $j - 1 ];
	}

	if ($a & $from_up) {
	    die "internal error" unless ($i > 0);
	    push @alt, [ $i - 1, $j ];
	}

	if ($a & $from_left) {
	    die "internal error" unless ($j > 0);
	    push @alt, [ $i, $j - 1];
	}

	if (!@alt) {
	    die "internal error";
	}

	my $cur = [ $i, $j ];
	my $move;
	if (@alt == 1) {
	    $move = $self->_simple_trace_back($cur, $alt[0],
					      $self->{callbacks});
	} else {
	    $move = $self->_trace_back($cur, \@alt);
	}

	if ($move eq 'align') {
	    --$i;
	    --$j;
	} elsif ($move eq 'shift_a') {
	    --$j;
	} elsif ($move eq 'shift_b') {
	    --$i;
	} else {
	    die "internal error";
	}
    }

    return $D->[$m]->[$n];
}

sub _align_affine {
    my $self = shift;
    my $a = shift;
    my $b = shift;

    my @D = ([ [ 0 ] ], [ [ 0 ] ], [ [ 0 ] ]); # indexed by $from_*_idx
    my $m = scalar(@$b);
    my $n = scalar(@$a);

    my $score_diag = sub {
        my ($i, $j) = @_;

	my @base = map { $_->[$i - 1]->[$j - 1]; } @D;
	my $base = max(@base);
	$base + &{$self->{score_sub}}($a->[$j - 1], $b->[$i - 1]);
    };

    my $score_up = sub {
	my ($i, $j) = @_;

	my @base = map { $_->[$i - 1]->[$j]; } @D;
	$base[$from_diag_idx] += $self->{gap_open_penalty};
	$base[$from_up_idx] += $self->{gap_extend_penalty};
	$base[$from_left_idx] += $self->{gap_open_penalty};
	max(@base);
    };

    my $score_left;
    if (!$self->{local}) {
	$score_left = sub {
	    my ($i, $j) = @_;

	    my @base = map { $_->[$i]->[$j - 1]; } @D;
	    $base[$from_diag_idx] += $self->{gap_open_penalty};
	    $base[$from_up_idx] += $self->{gap_open_penalty};
	    $base[$from_left_idx] += $self->{gap_extend_penalty};

	    max(@base);
	};
    } else {
	$score_left = sub {
	    my ($i, $j) = @_;

	    my @base = map { $_->[$i]->[$j - 1]; } @D;
	    if ($i < $m) {
		$base[$from_diag_idx] += $self->{gap_open_penalty};
		$base[$from_up_idx] += $self->{gap_open_penalty};
		$base[$from_left_idx] += $self->{gap_extend_penalty};
	    }

	    max(@base);
	};
    }

lib/Algorithm/NeedlemanWunsch.pm  view on Meta::CPAN

	foreach (@D) {
	    $_->[$i]->[0] = $self->{gap_open_penalty} +
		($i - 1) * $self->{gap_extend_penalty};
	}

	++$i;
    }

    # order must correspond to $from_* constants
    my @subproblems = ( $score_diag, $score_up, $score_left );

    $i = 1;
    while ($i <= $m) {
	$j = 1;
	while ($j <= $n) {
	    my $k = 0;
	    while ($k < 3) { # scalar(@D), scalar(@subproblems)
		$D[$k]->[$i]->[$j] = &{$subproblems[$k]}($i, $j);
		++$k;
	    }

	    # my $x = join ', ', map { $_->[$i]->[$j]; } @D;
	    # warn "$i, $j: $x\n";

	    ++$j;
	}

	++$i;
    }

    # like $score_up
    my @delta_up = ( $self->{gap_open_penalty}, $self->{gap_extend_penalty},
		     $self->{gap_open_penalty} );

    # like $score_left
    my @delta_left = ( $self->{gap_open_penalty}, $self->{gap_open_penalty},
		       $self->{gap_extend_penalty} );

    my @no_delta = (0, 0, 0);

    my @score = map { $_->[$m]->[$n]; } @D;
    my $res = max(@score);

    my $arrow = 0;
    my $flag = 1;
    my $idx = 0;
    while ($idx < 3) { # scalar(@score)
	if ($score[$idx] == $res) {
	    $arrow |= $flag;
	}

	$flag *= 2;
	++$idx;
    }

    $i = $m;
    $j = $n;
    while (($i > 0) || ($j > 0)) {
	my @alt;
	if ($arrow & $from_diag) {
	    die "internal error" unless ($i > 0) && ($j > 0);
	    push @alt, [ $i - 1, $j - 1 ];
	}

	if ($arrow & $from_up) {
	    die "internal error" unless ($i > 0);
	    push @alt, [ $i - 1, $j ];
	}

	if ($arrow & $from_left) {
	    die "internal error" unless ($j > 0);
	    push @alt, [ $i, $j - 1];
	}

	if (!@alt) {
	    die "internal error";
	}

	# my $x = join ', ', map { "[ " . $_->[0] . ", " . $_->[1] . " ]"; } @alt;
	# warn "$i, $j: $x\n";

	my $cur = [ $i, $j ];
	my $move;
	if (@alt == 1) {
	    $move = $self->_simple_trace_back($cur, $alt[0],
					      $self->{callbacks});
	} else {
	    $move = $self->_trace_back($cur, \@alt);
	}

	if ($move eq 'align') {
	    --$i;
	    --$j;

	    @score = map { $_->[$i]->[$j]; } @D;
	    if ($i == 0) {
		$arrow = $from_left;
	    } elsif ($j == 0) {
	        $arrow = $from_up;
	    } else {
		my $d = max(@score);
		$arrow = 0;
		$flag = 1;
		$idx = 0;
		while ($idx < 3) { # scalar(@score)
		    if ($score[$idx] == $d) {
			$arrow |= $flag;
		    }

		    $flag *= 2;
		    ++$idx;
		}
	    }
	} elsif ($move eq 'shift_a') {
	    --$j;

	    my @base = map { $_->[$i]->[$j] } @D;
	    my $delta;
	    if ($self->{local} && ($i == $m)) {
		$delta = \@no_delta;
	    } else {
		$delta = \@delta_left;
	    }

	    $arrow = $self->_retread($score[$from_left_idx],  $i, $j,
                \@base, $delta);
	    @score = @base;
	} elsif ($move eq 'shift_b') {
	    --$i;

	    my @base = map { $_->[$i]->[$j] } @D;
	    $arrow = $self->_retread($score[$from_up_idx], $i, $j,
		\@base, \@delta_up);
	    @score = @base;
	} else {
	    die "internal error";
	}
    }

    return $res;
}

sub _retread {
    my ($self, $to_score, $i, $j, $base, $delta) = @_;

    if ($i == 0) {
	return $from_left;
    } elsif ($j == 0) {
	return $from_up;
    }

    my $a = 0;
    my $flag = 1;
    my $idx = 0;
    while ($idx < 3) {
	if ($base->[$idx] + $delta->[$idx] == $to_score) {
	    $a |= $flag;
	}

	$flag *= 2;
	++$idx;
    }

    return $a;
}

sub _trace_back {
    my ($self, $cur, $sources) = @_;

    my $arg = { };
    foreach my $next (@$sources) {
        my $m = $self->_simple_trace_back($cur, $next, { });
	if ($m eq 'align') {
	    $arg->{align} = [ $cur->[1] - 1, $cur->[0] - 1 ];
	} elsif ($m eq 'shift_a') {
	    $arg->{shift_a} = $cur->[1] - 1;
	} elsif ($m eq 'shift_b') {
	    $arg->{shift_b} = $cur->[0] - 1;
	} else {
	    die "internal error";
	}
    }

    my $move;
    my $cb = $self->{callbacks};
    if (exists($cb->{select_align})) {
        $move = &{$cb->{select_align}}($arg);
	if (!exists($arg->{$move})) {
	    die "select_align callback returned invalid selection $move.";
	}
    } else {
        my @cn = qw(align shift_a shift_b);
	foreach my $m (@cn) {
	    if (exists($arg->{$m})) {
	        $move = $m;
		last;
	    }
	}

	if (!$move) {
	    die "internal error";
	}

	if (exists($cb->{$move})) {
	    if ($move eq 'align') {
	        &{$cb->{align}}(@{$arg->{align}});
	    } else {
	        &{$cb->{$move}}($arg->{$move});
	    }
	}
    }

    return $move;
}

sub _simple_trace_back {
    my ($self, $cur, $next, $cb) = @_;

    if ($next->[0] == $cur->[0] - 1) {
        if ($next->[1] == $cur->[1] - 1) {
	    if (exists($cb->{align})) {
	        &{$cb->{align}}($next->[1], $next->[0]);
	    }

	    return 'align';
	} else {
	    if ($next->[1] != $cur->[1]) {
	        die "internal error";
	    }

	    if (exists($cb->{shift_b})) {
	        &{$cb->{shift_b}}($cur->[0] - 1);
	    }

	    return 'shift_b';
	}
    } else {
        if ($next->[0] != $cur->[0]) {
	    die "internal error";
	}

	if ($next->[1] != $cur->[1] - 1) {
	    die "internal error";
	}

	if (exists($cb->{shift_a})) {
	    &{$cb->{shift_a}}($cur->[1] - 1);
	}

	return 'shift_a';
    }
}

1;

__END__

=head1 NAME

Algorithm::NeedlemanWunsch - sequence alignment with configurable scoring

=head1 VERSION

Version 0.03

=cut

=head1 SYNOPSIS

    use Algorithm::NeedlemanWunsch;

    sub score_sub {
        if (!@_) {
	    return -2; # gap penalty
        }

	return ($_[0] eq $_[1]) ? 1 : -1;
    }

    my $matcher = Algorithm::NeedlemanWunsch->new(\&score_sub);
    my $score = $matcher->align(
               \@a,
               \@b,
               {   align     => \&on_align,
                   shift_a => \&on_shift_a,
                   shift_b => \&on_shift_b,
		   select_align => \&on_select_align
               });

=head1 DESCRIPTION

Sequence alignment is a way to find commonalities in two (or more)
similar sequences or strings of some items or characters. Standard
motivating example is the comparison of DNA sequences and their
functional and evolutionary similarities and differences, but the
problem has much wider applicability - for example finding the longest
common subsequence (that is, C<diff>) is a special case of sequence
alignment.

Conceptually, sequence alignment works by scoring all possible
alignments and choosing the alignment with maximal score. For example,
sequences C<a t c t> and C<t g a t> may be aligned