package Align;

# Copyright (C) 2007 Martin A. Hansen.

# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

# http://www.gnu.org/copyleft/gpl.html


# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


# yak yak yak


# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


use strict;
use Data::Dumper;
use Storable qw( dclone );
use IPC::Open2;
use Maasha::Calc;
use Maasha::Seq;
use vars qw ( @ISA @EXPORT );

use constant {
    Q_BEG => 0,
    Q_END => 1,
    S_BEG => 2,
    S_END => 3,
    LEN   => 4,
    SCORE => 5,
    HEAD  => 0,
    SEQ   => 1,
};

@ISA = qw( Exporter );


# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


sub align_two_seq
{
    # Martin A. Hansen, August 2006.

    # Generates an alignment by chaining matches, which are subsequences
    # shared between two sequences. The routine functions by considering 
    # only matches within a given search space. If no matches are given
    # these will be generated, if long matches are found these will be
    # included in the alignment, otherwise matches will be included depending
    # on a calculated score. New search spaces spanning the spaces between
    # matches and the search space boundaries will be cast and recursed into.
    
    my ( $q_seq,     # sequence 1 ref
         $s_seq,     # sequence 2 ref
         $matches,   # list of matches
         $q_min,     # q sequence start position
         $q_max,     # q sequecne stop  position
         $s_min,     # s sequence start position
         $s_max,     # s sequecne stop  position
         $args,      # argument hash
       ) = @_;

    # returns a chain of matches that can be chained into an alignment

    $matches ||= [];
    $q_min   ||= 0;
    $s_min   ||= 0;
    $q_max   ||= length( ${ $q_seq } ) - 1;
    $s_max   ||= length( ${ $s_seq } ) - 1;

    $args->{ "long_matches" }  ||= 50;
    $args->{ "alph_len" }      ||= 4;

    my ( $wordsize, @chain, $match, $best_match, @long_matches );

    $matches = select_matches( $matches, $q_min, $q_max, $s_min, $s_max );

    if ( scalar @{ $matches } == 0 )   # no matches - find some!
    {
        # $wordsize = find_wordsize( $q_min, $q_max, $s_min, $s_max, $args );
        $wordsize = 4;
        $matches  = find_matches( $q_seq, $s_seq, $wordsize, $q_min, $q_max, $s_min, $s_max );

        while ( scalar @{ $matches } == 0 and $wordsize > 1 )
        {
            $wordsize--;
            $matches = find_matches( $q_seq, $s_seq, $wordsize, $q_min, $q_max, $s_min, $s_max );
        }
 
        if ( scalar @{ $matches } > 0 ) {
            push @chain, align_two_seq( $q_seq, $s_seq, $matches, $q_min, $q_max, $s_min, $s_max, $args );
        }
    }
    elsif ( @long_matches = grep { $_->[ LEN ] >= $args->{ "long_matches" } } @{ $matches } )   # long matches found - include all that don't overlap!
    {
        @long_matches = order_matches( \@long_matches );
        
        foreach $match ( @long_matches )
        {
            push @chain, $match;
            
            if ( $match->[ Q_BEG ] - $q_min >= 2 and $match->[ S_BEG ] - $s_min >= 2 ) { 
                push @chain, align_two_seq( $q_seq, $s_seq, $matches, $q_min, $match->[ Q_BEG ] - 1, $s_min, $match->[ S_BEG ] - 1, $args );   # intermediate search space
            }

            $q_min = $match->[ Q_END ] + 1;
            $s_min = $match->[ S_END ] + 1;
        }

        if ( $q_min + 1 < $q_max and $s_min + 1 < $s_max ) {
            push @chain, align_two_seq( $q_seq, $s_seq, $matches, $q_min, $q_max, $s_min, $s_max, $args );   # remaining search space
        }
    }
    else   # shorter matches are included according to score
    {
        foreach $match ( @{ $matches } ) {
            # $match->[ SCORE ] = score_match( $match, $q_min, $q_max, $s_min, $s_max );
            $match->[ SCORE ] = score_match_niels( $match, $q_seq, $q_min, $q_max, $s_min, $s_max );
        }

        # @{ $matches } = grep { $_->[ SCORE ] > 0 } @{ $matches };
        @{ $matches } = grep { $_->[ SCORE ] <= 0.25 } @{ $matches };
        # @{ $matches } = sort { $b->[ SCORE ] <=> $a->[ SCORE ] } @{ $matches };
        @{ $matches } = sort { $a->[ SCORE ] <=> $b->[ SCORE ] } @{ $matches };

        $best_match   = shift @{ $matches };

        if ( $best_match )
        {
            push @chain, $best_match;

            if ( $best_match->[ Q_BEG ] - $q_min >= 2 and $best_match->[ S_BEG ] - $s_min >= 2 ) {
                push @chain, align_two_seq( $q_seq, $s_seq, $matches, $q_min, $best_match->[ Q_BEG ] - 1, $s_min, $best_match->[ S_BEG ] - 1, $args );   # left search space
            }

            if ( $q_max - $best_match->[ Q_END ] >= 2 and $s_max - $best_match->[ S_END ] >= 2 ) {
                push @chain, align_two_seq( $q_seq, $s_seq, $matches, $best_match->[ Q_END ] + 1, $q_max, $best_match->[ S_END ] + 1, $s_max, $args );   # right search space
            }
        }
    }

    return wantarray ? @chain : \@chain;
}


sub select_matches
{
    # Martin A. Hansen, August 2006.

    # Given a list of matches and a search space,
    # include only those matches contained within
    # this search space.

    my ( $matches,   # list of matches
         $q_min,     # q sequence start position
         $q_max,     # q sequecne stop  position
         $s_min,     # s sequence start position
         $s_max,     # s sequecne stop  position
       ) = @_;

    # returns list of matches

    my ( @matches );

    @matches = grep { $_->[ Q_BEG ] >= $q_min and
                      $_->[ S_BEG ] >= $s_min and
                      $_->[ Q_END ] <= $q_max and
                      $_->[ S_END ] <= $s_max } @{ $matches };

    return wantarray ? @matches : \@matches;
}


sub order_matches
{
    # Martin A. Hansen, October 2006

    # given a list of long matches, order these by length and position
    # and include only those long matches that does not cross.

    my ( $long_matches,   # list of matches
       ) = @_;

    # returns a list of matches

    my ( @matches, $match, $i );

    @{ $long_matches } = sort { $b->[ LEN ] <=> $a->[ LEN ] } @{ $long_matches };
        
    @matches = shift @{ $long_matches };

    foreach $match ( @{ $long_matches } )
    {
        if ( $match->[ Q_END ] < $matches[ 0 ]->[ Q_BEG ] and $match->[ S_END ] < $matches[ 0 ]->[ S_BEG ] )
        {
            unshift @matches, $match;
        }
        elsif ( $match->[ Q_BEG ] > $matches[ -1 ]->[ Q_END ] and $match->[ S_BEG ] > $matches[ -1 ]->[ S_END ] )
        {
            push @matches, $match; 
        }
        else
        {
            for ( $i = 1; $i < @matches; $i++ )
            {
                if ( $matches[ $i - 1 ]->[ Q_END ] < $match->[ Q_BEG ] and $match->[ Q_END ] < $matches[ $i ]->[ Q_BEG ] and
                     $matches[ $i - 1 ]->[ S_END ] < $match->[ S_BEG ] and $match->[ S_END ] < $matches[ $i ]->[ S_BEG ]
                )
                {
                    splice @matches, $i, 0, dclone $match;
                    last;
                }
            }
        }
    }

    return wantarray ? @matches : \@matches; 
}


sub find_wordsize
{
    # Martin A. Hansen, August 2006.

    # Given a search space calculates the wordsize for a word so a match
    # occurs only a few times. More matches may be needed at low similarity in
    # order to avoid starting with a wrong match.

    my ( $q_min,   # q sequence start position
         $q_max,   # q sequecne stop position
         $s_min,   # s sequence start position
         $s_max,   # s sequecne stop position
         $args,    # argument hash
       ) = @_;

    # returns integer

    my ( $q_dim, $s_dim, $dim_min, $wordsize );

    $q_dim = $q_max - $q_min + 1;
    $s_dim = $s_max - $s_min + 1;

    $dim_min = Maasha::Calc::min( $q_dim, $s_dim );

    $wordsize = 1;

    if ( $dim_min > 2000000 )                # optimized for DNA
    {
        $wordsize = $args->{ "long_matches" };
    }
    elsif ( $dim_min > 100000 )                # optimized for DNA
    {
        $wordsize = int( $args->{ "long_matches" } / 2 );
    }
    elsif ( $q_dim > 100 or $s_dim > 100 )  # optimized for DNA
    {
        while ( $args->{ "alph_len" } ** $wordsize <= $q_dim * $s_dim and $wordsize < $dim_min ) {
            $wordsize++;
        }
    }
    else
    {
        while ( $args->{ "alph_len" } ** $wordsize <= $dim_min and $wordsize < $dim_min ) {
            $wordsize++;
        }
    }

    return $wordsize;
}


sub find_matches
{
    # Martin A. Hansen, November 2006

    # given two sequences, find all maximum expanded matches between these

    my ( $q_seq,      # sequence 1
         $s_seq,      # sequence 2
         $wordsize,   # word size
         $q_min,      # q sequence start position
         $q_max,      # q sequecne stop  position
         $s_min,      # s sequence start position
         $s_max,      # s sequecne stop  position
       ) = @_;

    # returns list of matches

    my ( $q_beg, $q_word, %word_hash, $s_beg, $s_word, $match, @matches );

    if ( length ${ $s_seq } > length ${ $q_seq } )
    {
        for ( $q_beg = $q_min; $q_beg <= $q_max - $wordsize + 1; $q_beg++ )
        {
            $q_word = lc substr ${ $q_seq }, $q_beg, $wordsize;

            next if $q_word =~ /n/i; #   DNA/genome optimization

            push @{ $word_hash{ $q_word } }, $q_beg;
        }

        for ( $s_beg = $s_min; $s_beg <= $s_max - $wordsize + 1; $s_beg++ )
        {
            $s_word = lc substr ${ $s_seq }, $s_beg, $wordsize;

            if ( exists $word_hash{ $s_word } )
            {
                foreach $q_beg ( @{ $word_hash{ $s_word } } )
                {
                    $match = [ $q_beg, $q_beg + $wordsize - 1, $s_beg, $s_beg + $wordsize - 1 ];

                    if ( grep { $match->[ Q_BEG ] >= $_->[ Q_BEG ] and
                                $match->[ Q_END ] <= $_->[ Q_END ] and
                                $match->[ S_BEG ] >= $_->[ S_BEG ] and
                                $match->[ S_END ] <= $_->[ S_END ] } @matches )
                    {
                        next;   # match is redundant
                    }
                    else
                    {
                        $match = expand_match( $q_seq, $s_seq, $match, $q_max, $q_min, $s_max, $s_min );
                        $match->[ LEN ] = $match->[ Q_END ] - $match->[ Q_BEG ] + 1;

                        push @matches, $match;
                    }
                }
            }
        }
    }
    else
    {
        for ( $s_beg = $s_min; $s_beg <= $s_max - $wordsize + 1; $s_beg++ )
        {
            $s_word = lc substr ${ $s_seq }, $s_beg, $wordsize;

            next if $s_word =~ /n/i; #   DNA/genome optimization

            push @{ $word_hash{ $s_word } }, $s_beg;
        }

        for ( $q_beg = $q_min; $q_beg <= $q_max - $wordsize + 1; $q_beg++ )
        {
            $q_word = lc substr ${ $q_seq }, $q_beg, $wordsize;

            if ( exists $word_hash{ $q_word } )
            {
                foreach $s_beg ( @{ $word_hash{ $q_word } } )
                {
                    $match = [ $q_beg, $q_beg + $wordsize - 1, $s_beg, $s_beg + $wordsize - 1 ];

                    if ( grep { $match->[ Q_BEG ] >= $_->[ Q_BEG ] and
                                $match->[ Q_END ] <= $_->[ Q_END ] and
                                $match->[ S_BEG ] >= $_->[ S_BEG ] and
                                $match->[ S_END ] <= $_->[ S_END ] } @matches )
                    {
                        next;   # match is redundant
                    }
                    else
                    {
                        $match = expand_match( $q_seq, $s_seq, $match, $q_max, $q_min, $s_max, $s_min );
                        $match->[ LEN ] = $match->[ Q_END ] - $match->[ Q_BEG ] + 1;

                        push @matches, $match;
                    }
                }
            }
        }
    }

    return wantarray ? @matches : \@matches;
}


sub expand_match
{
    # Martin A. Hansen, August 2006.

    # Given two sequences and a match, expand the match maximally.
    # A match is defined like this: [ Q_BEG, Q_END, S_BEG, S_END ]

    my ( $q_seq,   # sequence 1 ref
         $s_seq,   # sequence 2 ref
         $match,   # sequence match
         $q_max,   # q sequecne stop position
         $q_min,   # q sequence start position
         $s_max,   # s sequecne stop position
         $s_min,   # s sequence start position
       ) = @_;

    # returns match

    my ( $q_pos, $s_pos );

    # expanding forward

    $q_pos = $match->[ Q_END ] + 1;
    $s_pos = $match->[ S_END ] + 1;

    while ( $q_pos <= $q_max and $s_pos <= $s_max and lc substr( ${ $q_seq }, $q_pos, 1 ) eq lc substr( ${ $s_seq }, $s_pos, 1 ) )
    {
        $q_pos++;
        $s_pos++;
    }

    $match->[ Q_END ] = $q_pos - 1;
    $match->[ S_END ] = $s_pos - 1;

    # expanding backwards

    $q_pos = $match->[ Q_BEG ] - 1;
    $s_pos = $match->[ S_BEG ] - 1;

    while ( $q_pos >= $q_min and $s_pos >= $s_min and lc substr( ${ $q_seq }, $q_pos, 1 ) eq lc substr( ${ $s_seq }, $s_pos, 1 ) )
    {
        $q_pos--;
        $s_pos--;
    }

    $match->[ Q_BEG ] = $q_pos + 1;
    $match->[ S_BEG ] = $s_pos + 1;

    return $match;
}


sub score_match
{
    # Martin A. Hansen, August 2006
    
    # given a match and a search space scores the match according to three criteria:
    
    # 1) length of match - the longer the better.
    # 2) distance to closest corner - the shorter the better.
    # 3) distance to closest narrow end of the search space - the shorter the better.
    
    # each of these scores are divided by search space dimentions, and the total score
    # is calculated: total = score_len - score_corner - score_narrow
    
    # the higher the score, the better the match.
    
    my ( $match,   # match
         $q_min,   # q sequence start position
         $q_max,   # q sequecne stop  position
         $s_min,   # s sequence start position
         $s_max,   # s sequecne stop  position
       ) = @_;
    
    # returns a positive number
    
    my ( $q_dim, $s_dim, $score_len, $score_corner, $score_narrow, $score_total, $beg_diag_dist, $end_diag_dist,
         $min_diag_dist, $score_diag, $beg_narrow_dist, $end_narrow_dist, $max_narrow_dist );

    # ----- 1) scoring according to match length

    $score_len = $match->[ LEN ] ** 3;

    # ---- 2) score according to distance away from diagonal

    $q_dim = $q_max - $q_min + 1;
    $s_dim = $s_max - $s_min + 1;

    if ( $q_dim >= $s_dim ) # s_dim is the narrow end
    {
        $beg_diag_dist = Maasha::Calc::dist_point2line( $match->[ Q_BEG ], $match->[ S_BEG ], $q_min, $s_min, $q_min + $s_dim, $s_min + $s_dim );
        $end_diag_dist = Maasha::Calc::dist_point2line( $match->[ Q_BEG ], $match->[ S_BEG ], $q_max - $s_dim, $s_max - $s_dim, $q_max, $s_max );
    }
    else
    {
        $beg_diag_dist = Maasha::Calc::dist_point2line( $match->[ Q_BEG ], $match->[ S_BEG ], $q_min, $s_min, $q_min + $q_dim, $s_min + $q_dim );
        $end_diag_dist = Maasha::Calc::dist_point2line( $match->[ Q_BEG ], $match->[ S_BEG ], $q_max - $q_dim, $s_max - $q_dim, $q_max, $s_max );
    }

    $min_diag_dist = Maasha::Calc::min( $beg_diag_dist, $end_diag_dist );

    $score_diag = 2 * $min_diag_dist ** 2;

    # ----- 3) scoring according to distance to the narrow end of the search space

    if ( $q_dim > $s_dim ) # s_dim is the narrow end
    {
        $beg_narrow_dist = $match->[ Q_BEG ] - $q_min;
        $end_narrow_dist = $q_max - $match->[ Q_BEG ];

        $max_narrow_dist = Maasha::Calc::max( $beg_narrow_dist, $end_narrow_dist );
    }
    elsif ( $q_dim < $s_dim )
    {
        $beg_narrow_dist = $match->[ S_BEG ] - $s_min;
        $end_narrow_dist = $s_max - $match->[ S_BEG ];

        $max_narrow_dist = Maasha::Calc::max( $beg_narrow_dist, $end_narrow_dist );
    }
    else
    {
        $max_narrow_dist = 0;
    }

    $score_narrow = $max_narrow_dist;

    $score_total = $score_len - $score_narrow - $score_diag;
  #  $score_total = -1 if 3 * $min_diag_dist > $match->[ LEN ];

    return $score_total;
}


sub score_match_niels
{
    # Niels Larsen, June 2004.
    
    # Creates a crude "heuristic" attempt of telling how likely it is that a 
    # given match occurs by chance in a given search space. If sequences are
    # given their composition is taken into account. The scoring punishes 
    # distance from diagonal(s) and distance from previous match(es). Scores
    # range from zero and up, and lower is better. 

    my ( $match,   # Match array
         $q_seq,   # Either q_seq or s_seq
         $q_min,   # Lower bound search area (query sequence)
         $q_max,   # Upper bound search area (query sequence)
         $s_min,   # Lower bound search area (subject sequence)
         $s_max,   # Upper bound search area (subject sequence)
         ) = @_;
    
    # Returns a positive number. 
    
    my ( $q_beg, $s_beg, $q_end, $s_end, $q_dim, $s_dim, $seq, $pos,
         $q_delta_beg, $s_delta_beg, $q_delta_end, $s_delta_end, $i,
         $delta_beg_max, $delta_end_max, $as, $gs, $ts, $cs, $pmatch,
         $score, $moves, $dist_beg, $dist_end, $seqlen, %chars, $delta,
         $delta_max );
    
    $q_beg = $match->[Q_BEG];
    $q_end = $match->[Q_END];
    $s_beg = $match->[S_BEG];
    $s_end = $match->[S_END];
    
    # >>>>>>>>>>>>>>>>>>>>>>> CRUDE INITIAL SCORE <<<<<<<<<<<<<<<<<<<<<< 

    # Get the probability of a match from the sequence composition (when
    # match is longer than 20 and sequence is given, otherwise use 0.25)
    # and raise that to the power of the length. 

    if ( $match->[LEN] > 20 and defined $q_seq )
    {
        $seq = substr ${ $q_seq }, $q_beg, $q_end-$q_beg+1;
        $seqlen = length $seq;

        $chars{"a"} = $chars{"g"} = $chars{"c"} = $chars{"t"} = 0;

        for ( $i = 0; $i < $seqlen; $i++ )
        {
            $chars{ substr $seq, $i, 1 }++;
        }

        $pmatch = ($chars{"a"}/$seqlen)**2 + ($chars{"c"}/$seqlen)**2
                + ($chars{"g"}/$seqlen)**2 + ($chars{"t"}/$seqlen)**2;
    }
    else {
        $pmatch = 0.25;
    }

    $score = $pmatch ** ( $q_end - $q_beg + 1 ); 

#     # Punish by difference in height and width of search space,
    
#     $q_dim = $q_max - $q_min + 1;
#     $s_dim = $s_max - $s_min + 1;
    
#     if ( $q_dim != $s_dim ) {
#         $score *= abs ( $q_dim - $s_dim ) ** 2; 
#     }
    
    return $score if $score > 0.25;

    # Punish by how far the match is to the closest corner of the search
    # space,

    $q_delta_beg = $q_beg - $q_min;
    $s_delta_beg = $s_beg - $s_min;
    
    $q_delta_end = $q_max - $q_end;
    $s_delta_end = $s_max - $s_end;
    
    if ( $q_delta_beg > $s_delta_beg ) {
        $delta_beg_max = $q_delta_beg;
    } else {
        $delta_beg_max = $s_delta_beg;
    }

    if ( $q_delta_end > $s_delta_end ) {
        $delta_end_max = $q_delta_end;
    } else {
        $delta_end_max = $s_delta_end;
    }

    if ( $delta_beg_max <= $delta_end_max ) {
        $score *= ($delta_beg_max+1) ** 2.0;
    } else {
        $score *= ($delta_end_max+1) ** 2.0;
    }

    return $score if $score > 0.25;

    # Add penalty if the match is towards the narrow end of the
    # search space,

    if ( ($q_max - $q_min) <= ($s_max - $s_min) )
    {
        if ( $q_delta_beg > $s_delta_beg )
        {
            $score *= 2 * ( $q_delta_beg - $s_delta_beg ) ** 3;
        }
        elsif ( $q_delta_end > $s_delta_end )
        {
            $score *= 2 * ( $q_delta_end - $s_delta_end ) ** 3;
        }
    }
    else
    {
        if ( $s_delta_beg > $q_delta_beg )
        {
            $score *= 2 * ( $s_delta_beg - $q_delta_beg ) ** 3;
        }
        elsif ( $s_delta_end > $q_delta_end )
        {
            $score *= 2 * ( $s_delta_end - $q_delta_end ) ** 3;
        }        
    }

    if ( $score < 0 ) {
        print STDERR "q_min, q_max, s_min, s_max: $q_min, $q_max, $s_min, $s_max\n";
        die qq (Score <= 0 -> $score);
    }

    return $score;
}


sub print_alignment
{
    # Martin A. Hansen, August 2006.

    # Routine to print an alignment in fasta format based
    # on a list of matches and two given sequences.

    my ( $matches,  # list of matches
         $q_head,   # query   sequence head
         $q_seq,    # query   sequence ref
         $s_head,   # subject sequence head
         $s_seq,    # subject sequence ref
       ) = @_;

    my ( $q_pos, $s_pos, $q_nomatch, $q_match, $s_nomatch, $match, $q_aseq, $s_aseq, $i );

    @{ $matches } = sort { $a->[ Q_BEG ] <=> $b->[ Q_BEG ] } @{ $matches };

    $q_pos = 0;
    $s_pos = 0;

    for ( $i = 0; $i < @{ $matches }; $i++ )
    {
        $match = $matches->[ $i ];
    
        $q_nomatch = $match->[ Q_BEG ] - $q_pos;
        $s_nomatch = $match->[ S_BEG ] - $s_pos;

        if ( $q_nomatch - $s_nomatch > 0 ) 
        {
            $s_aseq .= "-" x ( $q_nomatch - $s_nomatch );
            $s_aseq .= substr ${ $s_seq }, $s_pos, $s_nomatch + $match->[ LEN ];
            $q_aseq .= substr ${ $q_seq }, $q_pos, $q_nomatch + $match->[ LEN ];
        }
        elsif ( $s_nomatch - $q_nomatch > 0 )
        {
            $q_aseq .= "-" x ( $s_nomatch - $q_nomatch );
            $q_aseq .= substr ${ $q_seq }, $q_pos, $q_nomatch + $match->[ LEN ];
            $s_aseq .= substr ${ $s_seq }, $s_pos, $s_nomatch + $match->[ LEN ];
        }
        else
        {
            $q_aseq .= substr ${ $q_seq }, $q_pos, $q_nomatch + $match->[ LEN ];
            $s_aseq .= substr ${ $s_seq }, $s_pos, $s_nomatch + $match->[ LEN ];
        }
    
        $q_pos += $q_nomatch + $match->[ LEN ];    
        $s_pos += $s_nomatch + $match->[ LEN ];    
    }

    $match = $matches->[ -1 ] || [ 0, 0, 0, 0, 0 ]; 

    $q_nomatch = length( ${ $q_seq } ) - $match->[ Q_END ];
    $s_nomatch = length( ${ $s_seq } ) - $match->[ S_END ];

    if ( $q_nomatch - $s_nomatch > 0 ) 
    {
        $q_aseq .= substr ${ $q_seq }, $q_pos, $q_nomatch + $match->[ LEN ];
        $s_aseq .= substr ${ $s_seq }, $s_pos, $s_nomatch + $match->[ LEN ];
        $s_aseq .= "-" x ( $q_nomatch - $s_nomatch );
    }
    elsif ( $s_nomatch - $q_nomatch > 0 )
    {
        $q_aseq .= substr ${ $q_seq }, $q_pos, $q_nomatch + $match->[ LEN ];
        $s_aseq .= substr ${ $s_seq }, $s_pos, $s_nomatch + $match->[ LEN ];
        $q_aseq .= "-" x ( $s_nomatch - $q_nomatch );
    }
    else
    {
        $q_aseq .= substr ${ $q_seq }, $q_pos, $q_nomatch + $match->[ LEN ];
        $s_aseq .= substr ${ $s_seq }, $s_pos, $s_nomatch + $match->[ LEN ];
    }

    print ">$q_head\n$q_aseq\n>$s_head\n$s_aseq\n";
}


sub print_matches
{
    # Martin A. Hansen, February 2007.

    my ( $matches,  # list of matches
         $q_head,   # query   sequence head
         $q_seq,    # query   sequence ref
         $s_head,   # subject sequence head
         $s_seq,    # subject sequence ref
         $args,     # argument hash         - OPTIONAL
       ) = @_;

    $args->{ "wrap" } ||= 80;

    my ( $q_pos, $s_pos, $match, $q_nomatch, $q_match, $s_nomatch, $q_aseq, $s_aseq, $pins, $i, $q, $s, $q_ruler, $s_ruler, $entries );

    @{ $matches } = sort { $a->[ Q_BEG ] <=> $b->[ Q_BEG ] } @{ $matches };

    $q_pos = 0;
    $s_pos = 0;

    for ( $i = 0; $i < @{ $matches }; $i++ )
    {
        $match = $matches->[ $i ];

        $q_nomatch = $match->[ Q_BEG ] - $q_pos;
        $s_nomatch = $match->[ S_BEG ] - $s_pos;

        $q = $q_pos;
        $s = $s_pos;

        if ( $q_nomatch - $s_nomatch > 0 ) 
        {
            $q_aseq .= substr ${ $q_seq }, $q_pos, ( $q_nomatch - $s_nomatch );
            $s_aseq .= "-" x ( $q_nomatch - $s_nomatch );
            $pins   .= " " x ( $q_nomatch - $s_nomatch );
            $q += ( $q_nomatch - $s_nomatch );
        }
        elsif ( $s_nomatch - $q_nomatch > 0 )
        {
            $s_aseq .= substr ${ $s_seq }, $s_pos, ( $s_nomatch - $q_nomatch );
            $q_aseq .= "-" x ( $s_nomatch - $q_nomatch );
            $pins   .= " " x ( $s_nomatch - $q_nomatch );
            $s += ( $s_nomatch - $q_nomatch );
        }

        while ( $q < $q_pos + $q_nomatch and $s < $s_pos + $s_nomatch )
        {
            $q_aseq .= substr ${ $q_seq }, $q, 1;
            $s_aseq .= substr ${ $s_seq }, $s, 1;

            if ( substr( ${ $q_seq }, $q, 1 ) eq substr( ${ $s_seq }, $s, 1 ) )
            {
                $pins .= ":";
            } else {
                $pins .= " ";
            }

            $q++;
            $s++;
        }

        $q_aseq .= substr ${ $q_seq }, $match->[ Q_BEG ], $match->[ LEN ];
        $s_aseq .= substr ${ $s_seq }, $match->[ S_BEG ], $match->[ LEN ];
        $pins   .= "|" x $match->[ LEN ];

        $q_pos += $q_nomatch + $match->[ LEN ];
        $s_pos += $s_nomatch + $match->[ LEN ];
    }

    $q_nomatch = length( ${ $q_seq } ) - ( $match->[ Q_END ] || 0 );
    $s_nomatch = length( ${ $s_seq } ) - ( $match->[ S_END ] || 0 );

    $q = $q_pos;
    $s = $s_pos;

    while ( $q < $q_pos + $q_nomatch and $q < length ${ $q_seq } and $s < $s_pos + $s_nomatch and $s < length ${ $s_seq } )
    {
        $q_aseq .= substr ${ $q_seq }, $q, 1;
        $s_aseq .= substr ${ $s_seq }, $s, 1;

        if ( substr( ${ $q_seq }, $q, 1 ) eq substr( ${ $s_seq }, $s, 1 ) ) {
            $pins .= ":";
        } else {
            $pins .= " ";
        }

        $q++;
        $s++;
        $q_pos++;
        $s_pos++;
    }

    if ( $q_nomatch - $s_nomatch > 0 ) 
    {
        $q_aseq .= substr ${ $q_seq }, $q_pos, ( $q_nomatch - $s_nomatch );
        $s_aseq .= "-" x ( $q_nomatch - $s_nomatch );
        $pins   .= " " x ( $q_nomatch - $s_nomatch );
    }
    elsif ( $s_nomatch - $q_nomatch > 0 )
    {
        $s_aseq .= substr ${ $s_seq }, $s_pos, ( $s_nomatch - $q_nomatch );
        $q_aseq .= "-" x ( $s_nomatch - $q_nomatch );
        $pins   .= " " x ( $s_nomatch - $q_nomatch );
    }

    $q_ruler = make_ruler( $q_aseq );    
    $s_ruler = make_ruler( $s_aseq );    

    $entries = [
        [ "ruler",     $q_ruler ],
        [ $q_head,     $q_aseq  ],
        [ "consensus", $pins    ],
        [ $s_head,     $s_aseq  ],
        [ "ruler",     $s_ruler ],
    ];

    align_print_multi( $entries, undef, $args->{ "wrap" } )
}


sub make_ruler
{
    # Martin A. Hansen, February 2007;

    # Gererates a ruler for a given sequence (with indels).

    my ( $seq
       ) = @_;

    # Returns string

    my ( $i, $char, $skip, $count, $gap, $tics );

    $gap = 0;
    $i   = 1;

    while ( $i <= length $seq )
    {
        $char = substr $seq, $i - 1, 1;
 
        $gap++ if $char eq "-";

        if ( $skip )
        {
            $skip--;
        }
        else
        {
            $count = $i - $gap;
            $count = 1 if $char eq "-";

            if ( $count % 100 == 0 )
            {
                if ( $count + length( $count ) >= length $seq )
                {
                    $tics .= "|";
                }
                else
                {
                    $tics .= "|" . $count;
                    $skip = length $count;
                }
            }
            elsif ( $count % 50 == 0 ) {
                $tics .= ":";
            } elsif ( $count % 10 == 0 ) {
                $tics .= ".";
            } else {
                $tics .= " ";
            }
        }

        $i++;
    }

    return $tics;
}


sub align_sim_local
{
    # Martin A. Hansen, May 2007.

    # Calculate the local similarity of an alignment based on
    # an alignment chain. The similarity is calculated as
    # the number of matching residues divided by the overall
    # length of the alignment chain. This means that a short
    # but "good" alignment will yield a high similarity, while
    # a long "poor" alignment will yeild a low similarity.

    my ( $chain,   # list of matches in alignment
       ) = @_;

    # returns a float

    my ( $match, $match_tot, $q_beg, $q_end, $s_beg, $s_end, $q_diff, $s_diff, $max, $sim );

    return 0 if not @{ $chain };

    $match_tot = 0;
    $q_end     = 0;
    $s_end     = 0;
    $q_beg     = 999999999;
    $s_beg     = 999999999;

    foreach $match ( @{ $chain } )
    {
        $match_tot += $match->[ LEN ];
    
        $q_beg = $match->[ Q_BEG ] if $match->[ Q_BEG ] < $q_beg;
        $s_beg = $match->[ S_BEG ] if $match->[ S_BEG ] < $s_beg;
    
        $q_end = $match->[ Q_END ] if $match->[ Q_END ] > $q_end;
        $s_end = $match->[ S_END ] if $match->[ S_END ] > $s_end;
    }

    $q_diff = $q_end - $q_beg + 1;
    $s_diff = $s_end - $s_beg + 1;

    $max = Maasha::Calc::max( $q_diff, $s_diff );

    $sim = sprintf( "%.2f", ( $match_tot / $max ) * 100 );

    return $sim;
}


sub align_sim_global
{
    # Martin A. Hansen, June 2007.

    # Calculate the global similarity of an alignment based on
    # an alignment chain. The similarity is calculated as
    # the number of matching residues divided by the
    # length of the shortest sequence.

    my ( $chain,   # list of matches in alignment
         $q_seq,   # ref to query sequence
         $s_seq,   # ref to subject sequence
       ) = @_;

    # returns a float

    my ( $match_tot, $min, $sim );

    return 0 if not @{ $chain };

    $match_tot = 0;

    map { $match_tot += $_->[ LEN ] } @{ $chain }; 

    $min = Maasha::Calc::min( length( ${ $q_seq } ), length( ${ $s_seq } ) );

    $sim = sprintf( "%.2f", ( $match_tot / $min ) * 100 );

    return $sim;
}


sub align_consensus
{
    # Martin A. Hansen, June 2006.

    # Given an alignment as a list of FASTA entries,
    # generates a consensus sequences based on the
    # entropies for each column similar to the way
    # a sequence logo i calculated. Returns the
    # consensus sequence as a FASTA entry.

    my ( $entries,   # list of aligned FASTA entries
         $type,      # residue type - OPTIONAL
         $min_sim,   # minimum similarity - OPTIONAL
       ) = @_;

    # Returns tuple

    my ( $bit_max, $data, $pos, $char, $score, $entry );

    $type    ||= Maasha::Seq::seq_guess_type( $entries->[ 0 ] );
    $min_sim ||= 50;

    if ( $type =~ /protein/ ) {
        $bit_max = 4;   
    } else {
        $bit_max = 2;
    }

    $data = Maasha::Seq::seqlogo_calc( $bit_max, $entries );

    foreach $pos ( @{ $data } )
    {
        ( $char, $score ) = @{ $pos->[ -1 ] };
        
        if ( ( $score / $bit_max ) * 100 >= $min_sim ) {
            $entry->[ SEQ ] .= $char;
        } else {
            $entry->[ SEQ ] .= "-";
        }
    }

    $entry->[ HEAD ] = "Consensus: $min_sim%";

    return wantarray ? @{ $entry } : $entry;
}


# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
