package Maasha::Match;

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


# Routines to match sequences


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


use strict;
use Data::Dumper;
use Storable qw( dclone );
use Maasha::Common;
use Maasha::Fasta;
use Maasha::Seq;
use Maasha::Berkeley_DB;
use vars qw ( @ISA @EXPORT );

@ISA = qw( Exporter );


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


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

    # Match sequences using MUMmer.

    my ( $entries1,   # FASTA entries
         $entries2,   # FASTA entries
         $options,    # additional MUMmer options - OPTIONAL
         $tmp_dir,    # temporary directory
       ) = @_;

    # Returns a list.

    my ( @args, $arg, $file_in1, $file_in2, $cmd, $file_out, $fh, $line, $result, @results );

    $tmp_dir ||= $ENV{ "BP_TMP" };

    $options->{ "word_size" } ||= 20;
    $options->{ "direction" } ||= "both";

    push @args, "-c";
    push @args, "-L";
    push @args, "-F";
    push @args, "-l $options->{ 'word_size' }";
    push @args, "-maxmatch";
    push @args, "-n" if not &Maasha::Seq::seq_guess_type( $entries1->[ 0 ]->[ 1 ] ) eq "protein";
    push @args, "-b" if $options->{ "direction" } =~ /^b/;
    push @args, "-r" if $options->{ "direction" } =~ /^r/;

    $arg = join " ", @args;

    $file_in1 = "$tmp_dir/muscle1.tmp";
    $file_in2 = "$tmp_dir/muscle2.tmp";
    $file_out = "$tmp_dir/muscle3.tmp";

    map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries1 };
    map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries2 };

    &Maasha::Fasta::put_entries( $entries1, $file_in1 );
    &Maasha::Fasta::put_entries( $entries2, $file_in2 );

    &Maasha::Common::run( "mummer", "$arg $file_in1 $file_in2 > $file_out 2>/dev/null" );

    $fh = &Maasha::Common::read_open( $file_out );

    while ( $line = <$fh> )
    {
        chomp $line;
        
        if ( $line =~ /^> (.+)Reverse\s+Len = (\d+)$/ )
        {
            $result->{ "Q_ID" }  = $1;
            $result->{ "Q_LEN" } = $2;
            $result->{ "DIR" }   = "reverse";
        }
        elsif ( $line =~ /^> (.+)Len = (\d+)$/ )
        {
            $result->{ "Q_ID" }  = $1;
            $result->{ "Q_LEN" } = $2;
            $result->{ "DIR" }   = "forward";
        }
        elsif ( $line =~ /^\s*(.\S+)\s+(\d+)\s+(\d+)\s+(\d+)$/ )
        {
            $result->{ "S_ID" }    = $1;
            $result->{ "S_BEG" }   = $2 - 1;
            $result->{ "Q_BEG" }   = $3 - 1;
            $result->{ "HIT_LEN" } = $4;
            $result->{ "S_END" }   = $result->{ "S_BEG" } + $result->{ "HIT_LEN" } - 1;
            $result->{ "Q_END" }   = $result->{ "Q_BEG" } + $result->{ "HIT_LEN" } - 1;

            push @results, dclone $result;
        }
    
    }

    unlink $file_in1;
    unlink $file_in2;
    unlink $file_out;

    return wantarray ? @results : \@results;
}


sub match_vmatch
{
    # Martin A. Hansen, April 2008.

    # Vmatches a list of records against a list of index files and the full
    # path to the result file is returned.

    my ( $tmp_dir,       # directory in where to save temp files
         $records,       # list of records
         $index_files,   # list of index files
         $options,       # argument hash
       ) = @_;

    # Returns a string.

    my ( $query_file, $result_file, @result_files, $fh_in, $fh_out, $line, @fields, $i, $record, $vmatch_args, @index_names, @seq_names, $count_list );

    $query_file  = "$tmp_dir/query.seq";
    $result_file = "$tmp_dir/vmatch.out";

    $fh_out = &Maasha::Common::write_open( $query_file );

    foreach $record ( @{ $records } ) 
    {
        if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
        {
            next if length $record->{ "SEQ" } < 12; # assuming that the index is created for 12 as minimum length

            push @seq_names, $record->{ "SEQ_NAME" };

            &Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ], $fh_out, 80 );
        }
    }

    close $fh_out;

    if ( $options->{ 'genome' } ) {
        $vmatch_args  = "-complete -d -p -q $query_file";
    } else {
        $vmatch_args  = "-complete -d -p -showdesc 100 -q $query_file";
    }

    $vmatch_args .= " -h " . $options->{ "hamming_dist" } if $options->{ "hamming_dist" };
    $vmatch_args .= " -e " . $options->{ "edit_dist" }    if $options->{ "edit_dist" };

    for ( $i = 0; $i < @{ $index_files }; $i++ )
    {
        &Maasha::Common::run( "vmatch", "$vmatch_args $index_files->[ $i ] > $result_file.$i" );

        push @result_files, "$result_file.$i";
    }

    unlink $query_file;

    $count_list = &vmatch_count_hits( \@result_files ) if ( $options->{ "count" } );

    $fh_out = &Maasha::Common::write_open( $result_file );

    for ( $i = 0; $i < @{ $index_files }; $i++ )
    {
        $index_files->[ $i ] =~ s/.+\/(.+)\.fna$/$1/ if $options->{ 'genome' };

        $fh_in = &Maasha::Common::read_open( "$result_file.$i" );
    
        while ( $line = <$fh_in> )
        {
            chomp $line;

            next if $line =~ /^#/;

            @fields = split " ", $line;

            next if $options->{ "max_hits" } and $count_list->[ $fields[ 5 ] ] > $options->{ 'max_hits' };

            $fields[ 1 ] = $index_files->[ $i ];                                     # S_ID
            $fields[ 9 ] = $count_list->[ $fields[ 5 ] ] if $options->{ "count" };   # SCORE
            $fields[ 5 ] = $seq_names[ $fields[ 5 ] ];                               # Q_ID

            print $fh_out join( "\t", @fields ), "\n";
        }

        close $fh_in;

        unlink "$result_file.$i";
    }

    close $fh_out;

    return $result_file;
}


sub vmatch_count_hits
{
    # Martin A. Hansen, April 2008.

    # Given a list of Vmatch result file, count duplications based
    # on q_id. The counts are returned in a list where the list index
    # corresponds to the q_id index in the query file.

    my ( $files,   # vmatch result files
       ) = @_;

    # Returns a list.

    my ( $file, $fh_in, $line, @fields, @count_list );

    foreach $file ( @{ $files } )
    {
        $fh_in = &Maasha::Common::read_open( $file );

        while ( $line = <$fh_in> )
        {
            chomp $line;

            next if $line =~ /^#/;

            @fields = split " ", $line;

            $count_list[ $fields[ 5 ] ]++;
        }

        close $fh_in;
    }

    return wantarray ? @count_list : \@count_list;
}


sub vmatch_count_hits_old
{
    # Martin A. Hansen, April 2008.

    # Given a Vmatch result file, substitute the
    # score field with the times the query sequence
    # was found.

    my ( $tmp_dir,     # directory in where to save temp files
         $path,        # full path to vmatch file
         $max_count,   # filter too abundant seqs - OPTIONAL
       ) = @_;

    # Returns nothing.

    my ( $fh_in, $fh_out, $line, @fields, @count_list );

    $fh_in = &Maasha::Common::read_open( $path );

    while ( $line = <$fh_in> )
    {
        chomp $line;

        next if $line =~ /^#/;

        @fields = split " ", $line;

        $count_list[ $fields[ 5 ] ]++;
    }

    close $fh_in;

    $fh_in  = &Maasha::Common::read_open( $path );
    $fh_out = &Maasha::Common::write_open( "$tmp_dir/vmatch.count" );

    while ( $line = <$fh_in> )
    {
        chomp $line;

        next if $line =~ /^#/;

        @fields = split " ", $line;

        $fields[ 9 ] = $count_list[ $fields[ 5 ] ];

        if ( $max_count ) {
            print $fh_out join( "\t", @fields ), "\n" if $fields[ 9 ] <= $max_count;
        } else {
            print $fh_out join( "\t", @fields ), "\n";
        }
    }

    close $fh_in;
    close $fh_out;

    rename "$tmp_dir/vmatch.count", $path;
}


sub vmatch_count_hits_old
{
    # Martin A. Hansen, April 2008.

    # Given a Vmatch result file, substitute the
    # score field with the times the query sequence
    # was found.

    my ( $tmp_dir,     # directory in where to save temp files
         $path,        # full path to vmatch file
         $max_count,   # filter too abundant seqs - OPTIONAL
       ) = @_;

    # Returns nothing.

    my ( $fh_in, $fh_out, $line, @fields, %count_hash );

    if ( $max_count ) {
        %count_hash = ();
    } else {
        %count_hash = &Maasha::Berkeley_DB::db_init( "$tmp_dir/hash.bdb" );
    }

    $fh_in = &Maasha::Common::read_open( $path );

    while ( $line = <$fh_in> )
    {
        chomp $line;

        next if $line =~ /^#/;

        @fields = split " ", $line;

        $count_hash{ $fields[ 5 ] }++;
    }

    close $fh_in;

    $fh_in  = &Maasha::Common::read_open( $path );
    $fh_out = &Maasha::Common::write_open( "$tmp_dir/vmatch.count" );

    while ( $line = <$fh_in> )
    {
        chomp $line;

        next if $line =~ /^#/;

        @fields = split " ", $line;

        $fields[ 9 ] = $count_hash{ $fields[ 5 ] };

        if ( $max_count ) {
            print $fh_out join( "\t", @fields ), "\n" if $fields[ 9 ] <= $max_count;
        } else {
            print $fh_out join( "\t", @fields ), "\n";
        }
    }

    close $fh_in;
    close $fh_out;

    if ( not $max_count )
    {
        untie %count_hash;
        unlink "$tmp_dir/hash.bdb";
    }

    rename "$tmp_dir/vmatch.count", $path;
}


sub vmatch_get_entry
{
    # Martin A. Hansen, January 2008.

    # Parses vmatch output records.

    my ( $fh,   # file handle to vmatch result file.
       ) = @_;

    # Returns a hash.

    my ( $line, @fields, %record );

    while ( $line = <$fh> )
    {
        chomp $line;

        next if $line =~ /^#/;

        @fields = split "\t", $line;

        $record{ "REC_TYPE" } = "VMATCH";

        $record{ "S_LEN" }      = $fields[ 0 ];
        $record{ "S_ID" }       = $fields[ 1 ];
        $record{ "S_BEG" }      = $fields[ 2 ];

        if ( $fields[ 3 ] eq "D" ) {
            $record{ "STRAND" } = "+";
        } else {
            $record{ "STRAND" } = "-";
        }

        $record{ "Q_LEN" }      = $fields[ 4 ];
        $record{ "Q_ID" }       = $fields[ 5 ];
        $record{ "Q_BEG" }      = $fields[ 6 ];
        $record{ "MATCH_DIST" } = $fields[ 7 ];
        $record{ "E_VAL" }      = $fields[ 8 ];
        $record{ "SCORE" }      = $fields[ 9 ];
        $record{ "IDENT" }      = $fields[ 10 ];

        $record{ "Q_END" }      = $record{ "Q_BEG" } + $record{ "Q_LEN" } - 1;
        $record{ "S_END" }      = $record{ "S_BEG" } + $record{ "S_LEN" } - 1;

        return wantarray ? %record : \%record;
    }
}


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


__END__
