#!/usr/bin/perl
# 
# pickBestTaxonomy.pl -- given multiple taxonomic analyses for the same set
#                of sequence data clustering, choose the "best" one for each
#                sequence, where best is defined by identifying the "deepest"
#                taxonomic classification.
#        All input files should be in the same format: tab or comma separated
#            table format (suitable for spreadsheet or database input);

# *************** This code assumes Ashley's "current" files are preprocessed
#                 to remove extraneous data and bad data: some of Ashley's
#               data files have incomplete names:  c__ or o__ when no name
#               has been found. These must be removed by the preprocessor.
      
#
# Input format: the header field is optional; try to infer if present, or force
#               it to have quoted strings, or force it to start with "seqId",
#               or require that it use the string "kingdom" in it somewhere, 
#               or use a command line option?
# 
#  header   |seqId, database, Kingdom, Phylum, Class, Order, Family, Genus, Species  
#  data     |seqid, dbIdent, k__name, p__name, c__name, o__name, g__name, s__name
#              however, any of the "name" fields can be missing

# 08/05/13 rdb v. 0.2
#              Fixed the incorrect field truncation in top-down version
#              Cleaned up the "bestDB" field so that a db id only gets included
#                 once when a non-empty best file is the input.
#              
# 07/22/13 rdb Created program: will read one db at a time, build hash table
#              with sequenceid as key and array of field data as value.
#              Strictly speaking, don't need hash table, since initial data
#              is coming in same order; easy, however, to not require that.
# 07/21/13 rdb: created

# 07/22/13 v. 0.1 sent to Ashley on 7/23
#

#------------------------------------------------------------------------------------

use strict;
use warnings;

my $usageMsg = q( Usage:  pickBestTaxonomy.pl [-b] bestTax.csv intax1.csv intax2.csv ... 
        where
            -b means to use "bottom-up" definition of best (lowest classification)
               default is "top-down": lowest complete classification from Kingdom down.
            *.csv are comma-separated taxonomy assignments from a database
                1st line is optional spreadsheet/db field names
                following lines:
                   seqId,dbId,Kingdom,Phylum,Class,Order,Family,Genus,Species.
            bestTax.csv need not exist; if it does, it is the first input file
                and the output file
            
    Program will update the bestTax.csv file with one line for each sequence that
    represents the "deepest" taxonomy for that sequence from all the inputs.
    
    Output files:
        bestTax.csv

        pickBest.log - will record any inconsistencies that are discovered.
                        if pickBest.log exists, will append to it.
                        LOG NOT YET IMPLEMENTED
    v. 0.2 08/05/13
);
#-------------------- Compile time constants ---------------------------
my $MAXBACKUP = 4;    # maximum number of backups to keep of the "best" file
#--------------------- main --------------------------------------------
if ( @ARGV < 2 || $ARGV[ 0 ] eq "-h" )
{
  die "$usageMsg\n";
}

my %bestTaxonomy = ();   # key is seqId, value is an array of fields: dbId, tax[]
my $header;             # will be determined by first file that has a header

my $topDown = 1;        # If defined, Best is the one with deepest first gap
my %topDownDepth = ();  #   else it is one with deepest specified type  
                          

if ( $ARGV[ 0 ] eq "-b" )
{
    $topDown = 0;
    shift( @ARGV );
}

my $bestFile = $ARGV[ 0 ];
shift( @ARGV );

if ( -e $bestFile )    # start out with best file
{
    &readFile( $bestFile );
    # ------- the following tries to avoid overwriting useful files by accident
    #         It renames recursively old files by appending ~.
    #         Probably dumb to do more than 2 or 3 levels
    if ( keys( %bestTaxonomy ) > 0 ) # if data in old best file, backup file
    {
        my $backup = "$bestFile~";
        my $bkupCount = 0;
        while ( -e $backup && $bkupCount < $MAXBACKUP )
        {
            rename( $backup, "$backup~" );
            $backup = "$backup~";
        }
        rename( $bestFile, "$bestFile~" );
    }
}

while ( @ARGV )
{
    &readFile( $ARGV[ 0 ] );
    shift( @ARGV );
}

&writeFile( $bestFile );

#------------------------- readFile( $filename  -------------------------------
# open and read the file specified;
#
# if a header has not been defined and one is present in this file, it will
#   become used to set the $header variable in the main
# Read each data line; and update the %bestTaxonomy hash table as appropriate.
#
sub readFile # ( $filename )
{
    my ( $filename ) = @_;
    my $error;
    my @buDepth = ( 0, 0, 0, 0, 0, 0, 0, 0 );  # bottom-up depth distr for file
    my @tdDepth = ( 0, 0, 0, 0, 0, 0, 0, 0 );  # top-down depth distr for file
    open IN, $filename or $error = 1;
    if ( defined( $error ))
    {
        print STDERR "Unable to open the file $filename\n";
        return -1;       
    }
    # skip over comments
    my $line = <IN>;
    while ( $line && substr( $line, 0, 1 ) eq "#" )
    {
        chomp( $line );
        $line = <IN>;
    }
    # handle first line; check if it's a proper header line
    if ( $line && $line =~ /Kingdom/ )
    {
        chomp( $line );
        if ( !defined( $header )) # if this is first header, save it
        {
            $header = $line;
            print "Updating header: $header\n";
        }
        $line = <IN>;
    }
    #
    # read "real" data lines
    while ( $line )
    {
        chomp( $line );
        if ( substr( $line, 0, 1 ) ne "#" )
        {
            my @fields = split( ",", $line );
            my $firstGap = 9;   # Species in array pos 8, so 9 is largest gap pos
            # for each field get rid of leading/trailing spaces
            for ( my $i = 0; $i < @fields; $i++ )
            {
                $fields[ $i ] =~ s/^\s*//; # trim leading spaces
                $fields[ $i ] =~ s/\s*$//; # trim trailing spaces
                if ( length( $fields[ $i ] ) == 0 && $i < $firstGap )
                {
                    $firstGap = $i;
                }
            }
            $buDepth[ @fields - 2 ]++;
            
            #-- get best data as an array -- if we have one for this sequence           
            my $seqId = $fields[ 0 ];
            my $bestData = $bestTaxonomy{ $seqId };
            my @bestArray = ();
            if ( $bestData )
            {
                @bestArray = @{ $bestData };
            }
            
            #-- if we are doing topdown version, delete extra fields
            if ( $topDown )
            {
                if ( $firstGap == 0 )  #  can only happen if there is bad input
                {
                    @fields = ();
                }
                else    # remove trailing fields
                {
                    @fields = splice( @fields, 0, $firstGap );  
                }
            }
            if ( @fields > @bestArray )
            {
                if ( $fields[ 0 ] =~ /denovo10003/ )
                {
                    print "------ @fields\n";
                }
                $bestTaxonomy{ $seqId } = \@fields;
            }
            elsif ( @fields == @bestArray )
            {
                my $dbId = $fields[ 1 ];
                if ( ! ( $bestArray[ 1 ] =~ /$dbId/ ))
                {
                    $bestArray[ 1 ] .= "+" . $fields[ 1 ];
                    $bestTaxonomy{ $seqId } = \@bestArray;
                }
            }
            $tdDepth[ @fields - 2 ]++;
        }        
        $line = <IN>;
    }
    print STDERR "$filename BU depth distribution: @buDepth\n";
    print STDERR "$filename TD depth distribution: @tdDepth\n";
}
#---------------------- writeFile -----------------------------------------
#
# If a header is defined, write the header, then write the data from
# hash table.
#
sub writeFile # ( outputFileName )
{
    my ( $outFile ) = @_;
    
    open OUT, ">$outFile" or die "Unable to open output file: $outFile";
    if ( $header )
    {
        print OUT "$header\n";
    }
    for my $key ( sort keys %bestTaxonomy )
    {
        my @fields = @{ $bestTaxonomy{ $key }}; # seems overly complicated
        my $line = join( ',', @fields );
        print OUT "$line\n";
    }
}