package Maasha::EMBL;

# 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 parse EMBL records.


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


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


@ISA = qw( Exporter );


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


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

    # Given a filehandle to an embl file,
    # fetches the next embl entry, and returns
    # this as a string with multiple lines.

    my ( $fh,    # filehandle to embl file
       ) = @_;

    # returns string
    
    my ( $entry );

    $/ = "//\n";

    $entry = <$fh>;

    return $entry;
}


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

    # given an embl entry extracts the keys
    # given in an argument hash. Special care
    # is taken to parse the feature table if
    # requested.

    my ( $entry,   # embl entry
         $args,    # argument hash
       ) = @_;

    # returns data structure

    my ( @lines, $line, %hash, $ft, $seq, $key );

    @lines = split "\n", $entry;

    foreach $line ( @lines )
    {
        if ( exists $args->{ "keys" } )
        {
            if ( $line =~ /^(\w{2})\s+(.*)/ and exists $args->{ "keys" }->{ $1 } )
            {
                if ( exists $hash{ $1 } and $1 eq "FT" ) {
                    $hash{ $1 } .= "\n" . $2;
                } elsif ( exists $hash{ $1 } ) {
                    $hash{ $1 } .= " " . $2;
                } else {
                    $hash{ $1 } = $2;
                }
            }
            elsif ( $line =~ /^\s+(.*)\s+\d+$/ and exists $args->{ "keys" }->{ "SEQ" } )
            {
                $seq .= $1;
            }
        }
        else
        {
            if ( $line =~ /^(\w{2})\s+(.*)/ )
            {
                if ( exists $hash{ $1 } and $1 eq "FT" ) {
                    $hash{ $1 } .= "\n" . $2;
                } elsif ( exists $hash{ $1 } ) {
                    $hash{ $1 } .= " " . $2;
                } else {
                    $hash{ $1 } = $2;
                }
            }
            elsif ( $line =~ /^\s+(.*)\s+\d+$/ )
            {
                $seq .= $1;
            }
        }
    }

    if ( $seq )
    {
        $seq =~ tr/ //d; 
        $hash{ "SEQ" } = $seq;
    }

#    foreach $key ( keys %hash )
#    {
#        next if $key =~ /^(SEQ|SEQ_FT|FT)/;
#
#        if ( not $hash{ $key } =~ /$args->{ $key }/i ) {
#            return wantarray ? () : {} ;
#        }
#    }

    if ( exists $hash{ "FT" } )
    {
        $seq =~ tr/ //d; 
        $ft = parse_feature_table( $hash{ "FT" }, $seq, $args );
        $hash{ "FT" } = $ft;
    }

    return wantarray ? %hash : \%hash;
}


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

    # parses the feature table of a EMBL/GenBank/DDBJ entry.
    # parsing takes place in 4 steps. 1) the feature key is
    # located. 2) the locator is located taking into # consideration
    # that it may be split over multiple lines, which is dealt with
    # by counting the params that always are present in multiline
    # locators. 3) the locator is used to fetch the corresponding
    # sequence. 4) qualifier key/value pars are located again taking
    # into consideration multiline values, which are dealt with by
    # keeping track of the "-count (value-less qualifers are also
    # included). only feature keys and qualifers defined in the
    # argument hash are returned.

    my ( $ft,     # feature table
         $seq,    # entry sequnce
         $args,   # argument hash
       ) = @_;

    # returns data structure

    my ( @lines, $key_regex, $i, $p, $q, %key_hash, $key, $locator, %qual_hash, $qual_name, $qual_val, $subseq );

    @lines = split "\n", $ft;

    $key_regex = "[A-Za-z0-9_']+";   # this regex should match every possible feature key (gene, misc_feature, 5'UTR ...)

    $i = 0;

    while ( $lines[ $i ] )
    {
        if ( $lines[ $i ] =~ /^($key_regex)\s+(.+)/ ) 
        {
            $key     = $1;
            $locator = $2;

            undef %qual_hash;

            # ---- getting locator

            $p = 1;

            if ( not balance_params( $locator ) )
            {
                while ( not balance_params( $locator ) )
                {
                    $locator .= $lines[ $i + $p ];
                    $p++;
                }
            }

            push @{ $qual_hash{ "_locator" } }, $locator;

            # ---- getting subsequence

            $subseq = parse_locator( $locator, $seq );

            push @{ $qual_hash{ "_seq" } }, $subseq;

            # ----- getting qualifiers

            while ( defined( $lines[ $i + $p ] ) and not $lines[ $i + $p ] =~ /^$key_regex/ )
            {
                if ( $lines[ $i + $p ] =~ /^\// )
                {
                    if ( $lines[ $i + $p ] =~ /^\/([^=]+)=(.*)$/ )
                    {
                        $qual_name = $1;
                        $qual_val  = $2;
                    }
                    elsif ( $lines[ $i + $p ] =~ /^\/(.*)$/ )
                    {
                        $qual_name = $1;
                        $qual_val  = "";
                    }

                    # ----- getting qualifier value

                    $q = 1;

                    if ( not balance_quotes( $qual_val ) )
                    {
                        while ( not balance_quotes( $qual_val ) )
                        {
                            $qual_val .= " " . $lines[ $i + $p + $q ];
                            $q++; 
                        }
                    }
                    
                    $qual_val =~ s/^"(.*)"$/$1/;
                    $qual_val =~ tr/ //d if $qual_name =~ /translation/i;
                    
                    if ( exists $args->{ "quals" } ) {
                        push @{ $qual_hash{ $qual_name } }, $qual_val if exists $args->{ "quals" }->{ $qual_name };
                    } else {
                        push @{ $qual_hash{ $qual_name } }, $qual_val;
                    }
                }
                    
                $p += $q;
            }

            if ( scalar keys %qual_hash > 0 )
            {
                if ( exists $args->{ "feats" } ) { 
                    push @{ $key_hash{ $key } }, dclone \%qual_hash if exists $args->{ "feats" }->{ $key };
                } else {
                    push @{ $key_hash{ $key } }, dclone \%qual_hash;
                }
            }
        }

        $i += $p;
    }

    return wantarray ? %key_hash : \%key_hash;
}


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

    # uses recursion to parse a locator string from a feature
    # table and fetches the appropriate subsequence. the operators
    # join(), complement(), and order() are handled.
    # the locator string is broken into a comma separated lists, and
    # modified if the params donnot balance. otherwise the comma separated
    # list of ranges are stripped from operators, and the subsequence are
    # fetched and handled according to the operators.
    # SNP locators are also dealt with (single positions).

    my ( $locator,   # locator string
         $seq,       # nucleotide sequence
         $subseq,    # subsequence being joined
         $join,      # join sequences
         $comp,      # complement sequence or not
         $order,     # order sequences
       ) = @_;

    # returns string

    my ( @intervals, $interval, $beg, $end, $newseq );

    @intervals = split ",", $locator;

    if ( not balance_params( $intervals[ 0 ] ) )                          # locator includes a join/comp/order of several ranges
    {
        if ( $locator =~ /^join\((.*)\)$/ )
        {
            $join = 1;
            $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
        }
        elsif ( $locator =~ /^complement\((.*)\)$/ )
        {
            $comp = 1;
            $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
        
        }
        elsif ( $locator =~ /^order\((.*)\)$/ )
        {
            $order = 1;
            $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
        }
    }
    else
    {
        foreach $interval ( @intervals )
        {
            if ( $interval =~ /^join\((.*)\)$/ )
            {
                $join = 1;
                $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
            }
            elsif ( $interval =~ /^complement\((.*)\)$/ )
            {
                $comp = 1;
                $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
            
            }
            elsif ( $interval =~ /^order\((.*)\)$/ )
            {
                $order = 1;
                $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
            }
            elsif ( $interval =~ /^[<>]?(\d+)[^\d]+(\d+)$/ )
            {
                $beg = $1;
                $end = $2;

                $newseq = substr $seq, $beg - 1, $end - $beg + 1;
    
                $newseq = Maasha::Seq::dna_revcomp( $newseq ) if $comp;

                if ( $order ) {
                    $subseq .= " " . $newseq;
                } else {
                    $subseq .= $newseq;
                }
            }
            elsif ( $interval =~ /^(\d+)$/ )
            {
                $beg = $1;

                $newseq = substr $seq, $beg - 1, 1 ;
    
                $newseq = Maasha::Seq::dna_revcomp( $newseq ) if $comp;

                if ( $order ) {
                    $subseq .= " " . $newseq;
                } else {
                    $subseq .= $newseq;
                }
            }
            else
            {
                warn qq(WARNING: Could not match locator -> $locator\n);
                # die qq(ERROR: Could not match locator -> $locator\n);
                $subseq .= "";
            }
        }
    }

    return $subseq;
}


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

    # given a string checks if left and right params
    # balances. returns 1 if balanced, else 0.

    my ( $string,   # string to check
       ) = @_;

    # returns boolean

    my ( $param_count );

    $param_count  = 0;
    $param_count += $string =~ tr/(//;
    $param_count -= $string =~ tr/)//;

    if ( $param_count == 0 ) {
        return 1;
    } else {
        return 0;
    }
}


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

    # given a string checks if the number of double quotes
    # balances. returns 1 if balanced, else 0.

    my ( $string,   # string to check
       ) = @_;

    # returns boolean

    my ( $quote_count );

    $quote_count = $string =~ tr/"//;

    if ( $quote_count % 2 == 0 ) {
        return 1;
    } else {
        return 0;
    }
}


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