package Maasha::AlignTwoSeq;

# 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 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,
    ALPH_LEN => 4,
};

@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
       ) = @_;

    # 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;

    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 );
        $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 );
        }
    }
    else   # matches are included according to score
    {
        foreach $match ( @{ $matches } ) {
            $match->[ SCORE ] = score_match( $match, $q_min, $q_max, $s_min, $s_max );
        }

        @{ $matches } = grep { $_->[ SCORE ] > 0 } @{ $matches };
        @{ $matches } = sort { $b->[ SCORE ] <=> $a->[ 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 );   # 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 );   # 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 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
       ) = @_;

    # 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 = 50;
    }
    elsif ( $dim_min > 100000 )             # optimized for DNA
    {
        $wordsize = 25;
    }
    elsif ( $q_dim > 100 or $s_dim > 100 )  # optimized for DNA
    {
        while ( ALPH_LEN ** $wordsize <= $q_dim * $s_dim and $wordsize < $dim_min ) {
            $wordsize++;
        }
    }
    else
    {
        while ( 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;

    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 insert_indels
{
    # Martin A. Hansen, June 2009.

    # Given a list of matches and references to q_seq and s_seq,
    # insert indels to form aligned sequences.

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

    # Returns nothing.

    my ( $q_indels, $s_indels, $match, $diff, $first );

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

    $first    = 1;
    $q_indels = 0;
    $s_indels = 0;

    # ---- initial indels ----

    $match = shift @{ $matches };

    substr ${ $q_seq }, $match->[ Q_BEG ] + $q_indels, $match->[ LEN ], uc substr ${ $q_seq }, $match->[ Q_BEG ] + $q_indels, $match->[ LEN ];
    substr ${ $s_seq }, $match->[ S_BEG ] + $s_indels, $match->[ LEN ], uc substr ${ $s_seq }, $match->[ S_BEG ] + $s_indels, $match->[ LEN ];

    $diff = ( $match->[ Q_BEG ] + $q_indels ) - ( $match->[ S_BEG ] + $s_indels );

    if ( $diff < 0 )
    {
        substr ${ $q_seq }, 0, 0, '-' x abs $diff;

        $q_indels += abs $diff;
    }
    elsif ( $diff > 0 )
    {
        substr ${ $s_seq }, 0, 0, '-' x abs $diff;

        $s_indels += abs $diff;
    }

    # ---- internal indels ----

    foreach $match ( @{ $matches } )
    {
        substr ${ $q_seq }, $match->[ Q_BEG ] + $q_indels, $match->[ LEN ], uc substr ${ $q_seq }, $match->[ Q_BEG ] + $q_indels, $match->[ LEN ];
        substr ${ $s_seq }, $match->[ S_BEG ] + $s_indels, $match->[ LEN ], uc substr ${ $s_seq }, $match->[ S_BEG ] + $s_indels, $match->[ LEN ];

        $diff = ( $match->[ Q_BEG ] + $q_indels ) - ( $match->[ S_BEG ] + $s_indels );

        if ( $diff < 0 )
        {
            substr ${ $q_seq }, $match->[ Q_BEG ] + $q_indels, 0, '-' x abs $diff;

            $q_indels += abs $diff;
        }
        elsif ( $diff > 0 )
        {
            substr ${ $s_seq }, $match->[ S_BEG ] + $s_indels, 0, '-' x abs $diff;

            $s_indels += abs $diff;
        }
    }

    # ---- last indels ----

    $diff = length( ${ $s_seq } ) - length( ${ $q_seq } );

    if ( $diff > 0 ) {
         ${ $q_seq } .= '-' x abs $diff;
    } else {
         ${ $s_seq } .= '-' x abs $diff;
    }
}


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;
}


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