#!/usr/bin/env perl

=head1 DESCRIPTION

Converts a JCVI metagenomic pipeline annotation file into a GO ontology gene association 
file.  The JCVI input files are columnar data with the following format:

    0   gene identifier
    1   gene product name label
    2   gene product name
    3   gene product name attribution
    4   gene symbol label
    5   gene symbol
    6   gene symbol attribution
    7   GO terms label
    8   GO terms
    9   GO terms attribution
    10  EC numbers label
    11  EC numbers
    12  EC numbers attribution
    13  CAZY identifier label
    14  CAZY identifier
    15  CAZY identifier attribution
    16  taxon label
    17  taxon
    18  taxon attribution
    

The format specification for the GAF is:

    http://www.geneontology.org/GO.format.gaf-2_0.shtml

Rather than populating the full GAF file, this script was initially written to just write enough 
to successfully run the map2slim binary.

The map2slim script can then be run like this:

map2slim /usr/local/projects/dacc/jcvi_metagenomic_autoannotate/output/PGA/HMPSlim.obo /usr/local/projects/dacc/jcvi_metagenomic_autoannotate/output/PGA/gene_ontology_ext.obo /usr/local/projects/dacc/jcvi_metagenomic_autoannotate/output/PGA/SRS011090/SRS011090.annotation.gaf

=head1 USAGE

convert_jcvi_annottab_to_gaf.pl /path/to/some_annotation.tab /path/to/output.gaf

The second argument, the output file, is optional.  If omitted, output will be
written to STDOUT and you can redirect it wherever you like.

=cut

use strict;
use warnings;
use POSIX;

my $tab_in = shift || die "ERROR: pass me an annotation.tab file or see perldoc";
my $out_file = shift || undef;

open(my $tab_fh, $tab_in) || die "failed to read annotation tab file: $!";

my $out_fh;

if ( $out_file ) {
    open($out_fh, ">$out_file") || die "failed to create output file: $!";
} else {
    $out_fh = *STDOUT;
}

my $file_date = strftime("%Y-%m-%d\n", localtime( (stat $tab_fh)[9] ));
chomp $file_date;

my @output_print_order = qw( 
    db               
    db_object_id     
    db_object_symbol 
    qualifier        
    go_id            
    db_reference     
    evidence_code    
    with_or_from     
    aspect           
    db_object_name   
    db_object_synonym
    db_object_type   
    taxon            
    date             
    assigned_by      
    annotation_ext   
    gene_prod_form_id
);

## this is the required first line of the output
print $out_fh "!gaf-version: 2.0\n";

my $gene_count_with_no_go_terms_assigned = 0;

while (my $line = <$tab_fh>) {
    chomp $line;
    
    my @cols = split("\t", $line);
    
    my @go_terms     = split(' \|\| ', $cols[8]);
    my @go_term_attr = split(' \|\| ', $cols[9]);
    
    for (my $go_i = 0; $go_i < scalar(@go_terms); $go_i++ ) {
    
        ## set defaults here
        my %output = (
            db                  => 'HMPDACC',
            db_object_id        => $cols[0],
            db_object_symbol    => $cols[5] ? $cols[5] : $cols[0],
            qualifier           => undef,               ## optional, according to spec
            go_id               => $go_terms[$go_i],
            db_reference        => $go_term_attr[$go_i],
            evidence_code       => 'IEA',
            with_or_from        => undef,               ## optional, according to spec
            aspect              => undef,               ## P, F or C
            db_object_name      => undef,               ## optional, according to spec
            db_object_synonym   => undef,               ## optional, according to spec
            db_object_type      => 'gene_product',
            taxon               => undef,
            date                => $file_date,
            assigned_by         => 'HMPDACC',
            annotation_ext      => undef,               ## optional, according to spec
            gene_prod_form_id   => undef,               ## optional, according to spec
        );

        
        
        ## output this line    
        for ( my $out_i=0; $out_i < scalar(@output_print_order); $out_i++ ) {
            defined $output{ $output_print_order[$out_i] } ? print $out_fh $output{ $output_print_order[$out_i] } : print $out_fh '';
            $out_i == scalar(@output_print_order) - 1 ? print $out_fh "\n" : print $out_fh "\t";
        }
    }
    
    $gene_count_with_no_go_terms_assigned++ if scalar(@go_terms) < 1;
}

print STDERR "File generated.  Gene count with no GO terms assigned: $gene_count_with_no_go_terms_assigned\n";






















