# $Id: FASTAH.pm,v 1.3 2007/01/09 17:16:00 rh4 Exp $
#
# BioMart module for BioMart::Formatter::FASTAH
#
# You may distribute this module under the same terms as perl
# itself.

# POD documentation - main docs before the code.

=head1 NAME

BioMart::Formatter::FASTAH

=head1 SYNOPSIS

FASTA formatter specific for Homologues

The FASTA Formatter returns whitespace separated tabular data
for a BioMart query's ResultTable
This is a FASTA formatter for Compara_homology which is very specific 
for the compara_mart project

=head1 DESCRIPTION

When given a BioMart::ResultTable containing the results of 
a BioMart::Query the FASTA Formatter will return tabular output
with one line for each row of data in the ResultTable and single spaces
separating the individual entries in each row. The getDisplayNames
and getFooterText can be used to return appropiately formatted
headers and footers respectively

=head1 AUTHORS

=over

=item *
Damian Smedley
=item *
Benoit Ballester
=back

=head1 CONTACT

This module is part of the BioMart project
http://www.biomart.org

Questions can be posted to the mart-dev mailing list:
mart-dev@ebi.ac.uk

=head1 METHODS

=cut

package BioMart::Formatter::FASTAH;

use strict;
use warnings;

# Extends BioMart::FormatterI
use base qw(BioMart::FormatterI);

# from CodonTable.pm in BioPerl
use vars qw(@NAMES @TABLES $TRCOL $CODONS %IUPAC_DNA 
	    $CODONGAP $GAP);# removed $TERMINATOR

# from CodonTable.pm in BioPerl
BEGIN { 
    use constant CODONSIZE => 3;
    #$GAP = '-';
    $CODONGAP = $GAP x CODONSIZE;

    @NAMES =			#id
	(
	 'Standard',		#1
	 'Vertebrate Mitochondrial',#2
	 'Yeast Mitochondrial',# 3
	 'Mold, Protozoan, and CoelenterateMitochondrial and Mycoplasma/Spiroplasma',#4
	 'Invertebrate Mitochondrial',#5
	 'Ciliate, Dasycladacean and Hexamita Nuclear',# 6
	 '', '',
	 'Echinoderm Mitochondrial',#9
	 'Euplotid Nuclear',#10
	 '"Bacterial"',# 11
	 'Alternative Yeast Nuclear',# 12
	 'Ascidian Mitochondrial',# 13
	 'Flatworm Mitochondrial',# 14
	 'Blepharisma Nuclear',# 15
	 'Chlorophycean Mitochondrial',# 16
	 '', '',  '', '',
	 'Trematode Mitochondrial',# 21
	 'Scenedesmus obliquus Mitochondrial', #22
	 'Thraustochytrium Mitochondrial' #23
	 );

    @TABLES =
	qw(
	   FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSS**VVVVAAAADDEEGGGG
	   FFLLSSSSYY**CCWWTTTTPPPPHHQQRRRRIIMMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSSSVVVVAAAADDEEGGGG
	   FFLLSSSSYYQQCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   '' ''
	   FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNNKSSSSVVVVAAAADDEEGGGG
	   FFLLSSSSYY**CCCWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   FFLLSSSSYY**CC*WLLLSPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSGGVVVVAAAADDEEGGGG
	   FFLLSSSSYYY*CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNNKSSSSVVVVAAAADDEEGGGG
	   FFLLSSSSYY*QCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   FFLLSSSSYY*LCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   '' '' '' ''
	   FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNNKSSSSVVVVAAAADDEEGGGG   
	   FFLLSS*SYY*LCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   FF*LSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
	   );
  

    my @nucs = qw(t c a g);
    my $x = 0;
    ($CODONS, $TRCOL) = ({}, {});
    for my $i (@nucs) {
	for my $j (@nucs) {
	    for my $k (@nucs) {
		my $codon = "$i$j$k";
		$CODONS->{$codon} = $x;
		$TRCOL->{$x} = $codon;
		$x++;
	    }
	}
    }

    %IUPAC_DNA = ( A => [qw(A)],
         C => [qw(C)],
         G => [qw(G)],
         T => [qw(T)],
         U => [qw(U)],
         M => [qw(A C)],
         R => [qw(A G)],
         W => [qw(A T)],
         S => [qw(C G)],
         Y => [qw(C T)],
         K => [qw(G T)],
         V => [qw(A C G)],
         H => [qw(A C T)],
         D => [qw(A G T)],
         B => [qw(C G T)],
         X => [qw(G A T C)],
         N => [qw(G A T C)]
         );
    #$TERMINATOR = '*';
}
#--- END from bioperl / GenomicSequence.pm

sub _new {
    my ($self) = @_;

    $self->SUPER::_new();
}

sub processQuery {
    my ($self, $query) = @_;

    $self->set('original_attributes',[@{$query->getAllAttributes()}]) 
	if ($query->getAllAttributes());
    $self->set('query',$query);
    return $query;
}

sub nextRow {
    my $self = shift;

    my $rtable = $self->get('result_table');
    my $row = $rtable->nextRow;
    if (!$row){
        return;
    }
#    my $array_length = @{$row};
    
    #-- determine  which coordinate we are using
    my ($tag1, $tag2);
    my (@coorarray1, @coorarray2);

    my $k = 0;
    my $j = 0;
    foreach my $attribute  (@$row){
	#print "attribute: $attribute\n";
	if (not defined $attribute){$attribute = "NULL"}
	if ($attribute =~ m/tagExonseq|tagCdna|tagCodingseq|tagUnsp|tag3utr|tag5utr|tagPeptide/){
	    $j++;
	    if ($j == 1){
		#print " attribute1::  $attribute\n";
		@coorarray1 = split(/\|/, $attribute);
		$tag1  = shift (@coorarray1); # remove the tags
	    }
	    if ($j == 2){
		#print " attribute2::  $attribute\n";
		@coorarray2 = split(/\|/, $attribute);
		$tag2  = shift (@coorarray2);
	    }
	  }  
	$k++;
    } 
#    print "TAGS : $tag1 $tag2\n";
    if ($tag1 ne $tag2){return "Tag are not the same:".$tag1." - ".$tag2."\n\n";}
    
    #-- then format depending of the tag (cdna, peptide, etc.)
    if ($tag1 =~ m/tagExonseq|tagCdna|tagCodingseq|tagUnsp|tag3utr|tag5utr/){

	&formatseq($row);

    }elsif ($tag1  =~ m/tagPeptide/){
	
	&formatseqPep($row);
#	print "tagPeptide : $tag1 $tag2\n";
    }else {
	return "tags not supported by the formater FASTAH.pm  tags(tag1:$tag1 tag2:$tag2)\n";
   }

}
#------------------
sub formatseq {
    my ($row)= shift;
    my ($finalseq1, $finalseq2, $header_atts1, $header_atts2, 
	$strand1, $strand2, $tag1, $tag2) ;
    my @atts1 = ();
    my @atts2 = ();
    my @coorarray = ();
    my @coorarray1;
    my @coorarray2;

    
    #--  get sequences
    my $seq1 = shift @$row;
    my $seq2 = shift @$row;

#    my $i = 0 ;
#    foreach my $attribute  (@$row){
#	print "attribute: $attribute\n";
#	if (!$attribute){$attribute = "NULL"}
#	elsif ($attribute eq "tagHeader"){
#	    push (@atts1, $row->[$i+1]);
#	    push (@atts2, $row->[$i+2]);
#	}
#	$i++;
#    }


    my $i = 0 ;
    my $k = 0;
    my $j = 0;
    foreach my $attribute  (@$row){
	#print "attribute: $attribute\n";
	#-- Get the attributes 
	if (!$attribute){$attribute = "NULL"}
	if ($attribute eq "tagHeader"){
	    push (@atts1, $row->[$i+1]);
	    push (@atts2, $row->[$i+2]);
	}
	$i++;
	
	#-- Get the coordinates
	if ($attribute =~ /^(tagExonseq|tagCdna|tagCodingseq|tagUnsp|tag3utr|tag5utr)/){
	    $j++;
	    if ($j == 1){
#		print " attribute1::  $attribute\n";
		@coorarray1 = split(/\|/, $attribute);
		$strand1 = $row->[$k+2];
	    }
	    if ($j == 2){
#		print " attribute2::  $attribute\n";
		@coorarray2 = split(/\|/, $attribute);
		$tag1  = shift (@coorarray1); # remove the tags
		$tag2  = shift (@coorarray2);
		$strand2 = $row->[$k+2];
	    }
	  }  
	$k++;
    } 
    
    #-- Chop the seq according to the given coordinate
    foreach my $coors1 (@coorarray1){
	if ($coors1 eq "No UTR"){
	    $finalseq1 = "No UTR";
	}else{
	    my @coor1 = split(/:/,$coors1);
	    $finalseq1 .= substr ($seq1, $coor1[0], $coor1[1]);
	}
    }
    
    foreach my $coors2 (@coorarray2){
	if ($coors2 =~ m/No UTR/gi){
	    $finalseq2 = "No UTR";
	}else{
	    my @coor2 = split(/:/,$coors2);
	    $finalseq2 .= substr ($seq2, $coor2[0], $coor2[1]);
	}
    }
    #-- RC the seq if needed !
    if (($strand1 < 0)&&($finalseq1 ne "No UTR")) {$finalseq1 = &rc($finalseq1);}
    if (($strand2 < 0)&&($finalseq2 ne "No UTR")) {$finalseq2 = &rc($finalseq2);}
	
    #-- format the seq whith a 60bp width
    $finalseq1 =~ s/(\w{60})/$1\n/g;
    $finalseq2 =~ s/(\w{60})/$1\n/g;

    #-- Join the attributes together
    $header_atts1 = join ("|", @atts1);
    $header_atts2 = join ("|", @atts2);
    
    #--  Then return the whole thing
    return ">" . $header_atts1 . "\n"
	       . $finalseq1 ."\n"
	       . ">" . $header_atts2 . "\n"
	       . $finalseq2 ."\n"
	       . "#\n";

}
#------------------
sub formatseqPep {
    my ($row)= shift;
    my ($finalseq1, $finalseq2, $header_atts1, $header_atts2, 
	$strand1, $strand2, $tag1, $tag2, $chro1, $chro2,
	$seqedits1, $seqedits2) ;
    my @atts1 = ();
    my @atts2 = ();
    my @coorarray = ();
    my @coorarray1;
    my @coorarray2;
    
    #--  get sequences
    my $seq1 = shift @$row;
    my $seq2 = shift @$row;

    my $i = 0 ;
    my $k = 0;
    my $j = 0;
    foreach my $attribute  (@$row){
	#print "attribute: $attribute\n";
	#-- Get the attributes 
	if (!$attribute){$attribute = "NULL"}
	if ($attribute eq "tagHeader"){
	    push (@atts1, $row->[$i+1]);
	    push (@atts2, $row->[$i+2]);
	}
	$i++;
	
	#-- Get the coordinates
	if ($attribute =~ /^(tagPeptide)/){
	    $j++;
	    if ($j == 1){
#		print " attribute1::  $attribute\n";
		@coorarray1 = split(/\|/, $attribute);
		$strand1 = $row->[$k+2];
		$seqedits1 = $row->[$k+4];
		$chro1 = $row->[$k+6];
		
	    }
	    if ($j == 2){
#		print " attribute2::  $attribute\n";
		@coorarray2 = split(/\|/, $attribute);
		$tag1  = shift (@coorarray1); # remove the tags
		$tag2  = shift (@coorarray2);
		$strand2 = $row->[$k+2];
		$seqedits2 = $row->[$k+4];
		$chro2 = $row->[$k+6];
	    }
	  }  
	$k++;
    } 

#    print "tag1:$tag1 strand1:$strand1 seqedits1:$seqedits1  chro1:$chro1\n";
#    print "tag2:$tag2 strand2:$strand2 seqedits2:$seqedits2  chro2:$chro2\n";

    
    #-- Chop the seq according to the given coordinate
    foreach my $coors1 (@coorarray1){
	    my @coor1 = split(/:/,$coors1);
	    $finalseq1 .= substr ($seq1, $coor1[0], $coor1[1]);
    }
    
    foreach my $coors2 (@coorarray2){
	    my @coor2 = split(/:/,$coors2);
	    $finalseq2 .= substr ($seq2, $coor2[0], $coor2[1]);
    }


    #-- RC the seq if needed !
    if ($strand1 < 0) {$finalseq1 = &rc($finalseq1);}
    if ($strand2 < 0) {$finalseq2 = &rc($finalseq2);}
	
    #-- Translated it
    $finalseq1 = &translate($finalseq1,$chro1);
    $finalseq2 = &translate($finalseq2,$chro2);

    #-- Post process the seq with seqedits
    if ($finalseq1 =~ /\*$/){ chop $finalseq1 }
    if ($finalseq2 =~ /\*$/){ chop $finalseq2 }
    
    #-- Apply seq edits if exists
    if ($seqedits1  ne 'NULL'){
	my @seqed = split (/\|/,$seqedits1);
	shift @seqed; #remove tag
	foreach my $se (@seqed) {
	#    $finalseq1 =  
	    &apply_edit($se, \$finalseq1);
	}
    }
     if ($seqedits2  ne 'NULL'){
	my @seqed = split (/\|/,$seqedits2);
	shift @seqed; #remove tag
	foreach my $se (@seqed) {
	    #$finalseq2 =  
	    &apply_edit($se, \$finalseq2);
	}
    }
    
    #-- format the seq whith a 60bp width
    $finalseq1 =~ s/(\w{60})/$1\n/g;
    $finalseq2 =~ s/(\w{60})/$1\n/g;

    #-- Join the attributes together
    $header_atts1 = join ("|", @atts1);
    $header_atts2 = join ("|", @atts2);
    
    #--  Then return the whole thing
    return ">" . $header_atts1 . "\n"
	       . $finalseq1 ."\n"
	       . ">" . $header_atts2 . "\n"
	       . $finalseq2 ."\n"
	       . "#\n";

}

#----------------------
sub rc{
    my ($seq) = @_;

    $seq = reverse($seq);
    $seq =~ tr/YABCDGHKMRSTUVyabcdghkmrstuv/RTVGHCDMKYSAABrtvghcdmkysaab/;

    return $seq;
}
#----------------------
#----------------------
sub getDisplayNames {
    my $self = shift;
    return '';
}

sub isSpecial {
    return 1;
}
#------------------
sub apply_edit  {
    my $seqeds   = shift;
    my $seqref = shift;
    my @seqedits = split (/:/,$seqeds);
    
    #if(ref($seqref) ne 'SCALAR') {
    #  throw("Reference to scalar argument expected");
    #}
    
    if(!defined($seqedits[0]) || !defined($seqedits[1])) {
	return $seqref;
    }
    
    my $len = $seqedits[1] - $seqedits[0] + 1;
    substr($$seqref, $seqedits[0] - 1, $len) = $seqedits[2];
    
    return $seqref;
}
#-----------------------------------------
#-------From Bioperl CodonTable.pm -------
#-------Modified like GenomicSequence ----
sub translate {
    my ($seq, $chro) = @_;
    #warn("Calling translate without a seq argument!") unless defined $seq;
    return '' unless $seq;

    #my $id = $self->id;
    # 1 = vertebrate  2= Vertebrate Mitochondrial
    my $id;#my $id = '1'; # by default it's1
    if ($chro =~ /MT/g){$id = '2'}else{$id = '1'}
 
    my ($partial) = 0;
    $partial = 2 if length($seq) % CODONSIZE == 2;
    
    $seq = lc $seq;
    $seq =~ tr/u/t/;
    my $protein = "";
    if ($seq =~ /[^actg]/ ) { #ambiguous chars
        for (my $i = 0; $i < (length($seq) - (CODONSIZE-1)); $i+= CODONSIZE) {
            my $triplet = substr($seq, $i, CODONSIZE);
#	    if( $triplet eq $CODONGAP ) {
#		$protein .= $GAP;
#	    } els
	    if (exists $CODONS->{$triplet}) {
		$protein .= substr($TABLES[$id-1], 
				   $CODONS->{$triplet},1);
	    } else {
		$protein .= _translate_ambiguous_codon($triplet, $id);
	    }
	}
    } else { # simple, strict translation
	for (my $i = 0; $i < (length($seq) - (CODONSIZE -1)); $i+=CODONSIZE) {
            my $triplet = substr($seq, $i, CODONSIZE); 
#            if( $triplet eq $CODONGAP ) {
#		$protein .= $GAP;
#	    } 
	    if (exists $CODONS->{$triplet}) {
                $protein .= substr($TABLES[$id-1], $CODONS->{$triplet}, 1);
	    } else {
                $protein .= 'X';
            }
        }
    }
    if ($partial == 2) { # 2 overhanging nucleotides
	my $triplet = substr($seq, ($partial -4)). "n";
#	if( $triplet eq $CODONGAP ) {
#	    $protein .= $GAP;
#	} elsif (exists $CODONS->{$triplet}) {
	
	if (exists $CODONS->{$triplet}) {
	    my $aa = substr($TABLES[$id-1], $CODONS->{$triplet},1);       
	    $protein .= $aa;
	} else {
	    $protein .= _translate_ambiguous_codon($triplet, $id, $partial);
	}
    }
    return $protein;
}
#----------------------
sub _translate_ambiguous_codon {
    my ($triplet, $id, $partial) = @_;
    $partial ||= 0;
    #my $id = '1';
    my $aa;
    my @codons = _unambiquous_codons($triplet);
    my %aas =();
    foreach my $codon (@codons) {
	$aas{substr($TABLES[$id-1],$CODONS->{$codon},1)} = 1;
    }
    my $count = scalar keys %aas;
    if ( $count == 1 ) {
	$aa = (keys %aas)[0];
    }
    elsif ( $count == 2 ) {
	if ($aas{'D'} and $aas{'N'}) {
	    $aa = 'B';
	}
	elsif ($aas{'E'} and $aas{'Q'}) {
	    $aa = 'Z';
	} else {
	    $partial ? ($aa = '') : ($aa = 'X');
	}
    } else {
	$partial ? ($aa = '') :  ($aa = 'X');
    }
    return $aa;
}
#----------------------
sub _unambiquous_codons{
    my ($value) = @_;
    my @nts = ();
    my @codons = ();
    my ($i, $j, $k);
    @nts = map { $IUPAC_DNA{uc $_} }  split(//, $value);
    for my $i (@{$nts[0]}) {
	for my $j (@{$nts[1]}) {
	    for my $k (@{$nts[2]}) {
		push @codons, lc "$i$j$k";
	    }
	}
    }
    return @codons;
}

1;
