#!/usr/bin/perl -w
use strict;
use Getopt::Long;

# annspec2gmm.pl
# 04 May 2004
# originally written: 17 March 2004, Gordon Robertson
#
# converts ANN-Spec output file into ;-delimited GMM-format files for 2D viewer
# retains input file order, and knows seq length without reading MFA
#
# 09avr04: agr: changed motif basename to 'ann_r' so motif names in HitPlotter are more meaningful
#
# GMM file format
# Seq_name;seq_length;block1;block2;...
# ...where each block describes a 'hit' with:
#    matrix_name;relstart;relend;score;strand
# ...and relstart and relend are genomic coordinates of a hit relative to the submitted FASTA sequence that had the hit
#    counting from the left end of the sequence, with all sequences being transcribed from left to right
#
# Example command line arguments:
# -ann=annspec_W8.results -out=test_annspec_W8.gmm
#
#
# Required arguments: 
my ( $results_in, $gmm_outfile ); # read in files
GetOptions( "ann=s" => \$results_in,
	    "out=s" => \$gmm_outfile
	    );
open IN,  "< $results_in" or die "Can't open $results_in !\n";
open OUT, "> $gmm_outfile" or die "Can't open $gmm_outfile !\n";	

my ( $scoring_rank, $training_run, $train_iter_best_score_found );
my ( @seq_names, $seq_name, $seq_length, $input_file            );
my ( @comparison_array, $tot_seqs ); # this is the GMM-formatted data
my ( $matrix_name, $best_score, $weight_cons, $alignment_cons, $gene_name_entry, $this_seq_hitdata );
my ( $n_site_index, $S_wt_mtrx_site_score, $Lpos_Lend_offset, $relend, $width, %hdata     );
my ( $complement_strand_flag, $sites, $strand, $Rpos_Rend_offset, $hitdata, $info_content );    
    
my $seqnbr     = 1;
my $matrix_nbr = 0; #this is actually rank number...

# Retain the MFA's sequence order in the GMM output file
# ...and so in the order of sequences in the 2D viewer
my $line;
READLINE: while ( my $line = <IN> ) {
    # Example:
    # section 1: gene names and indexes, seq lengths
    # sequence 1: F36D3.9
    #      Gene Index:     1
    #      Ortholog Index: 1
    #      fragments: 1-1000 1-1000
    # sequence 2: CBG23343
    #      Gene Index:     1
    #      Ortholog Index: 2
    #      fragments: 1-1000 1-1000
    # sequence 3: C15C8.3
    #      Gene Index:     2
    #      Ortholog Index: 1
    #      fragments: 1-865 1-865
    #...
    #
    # E.g. SQI "../input/adult_66pairs_1500bp_genes.mfa"
    if ( $line =~ /^SQI "([^"]+)/ ) {
	$input_file = $1;
    }

    # E.g.: SQI 1 >H22K11.1(asp-3)	(1-1500, A:522, C:210, G:246, T:522)
    if ( $line =~ /^SQI\s(\d+)\s>(\S+)\s\(1-([^\,]+)/ ) {
        $seqnbr =     $1;
        $seq_name =   $2;
	$seq_length = $3;
        $seq_names[$seqnbr -1] = $seq_name; # loading an array, so can keep sequence order later!
        print "sequence: $seqnbr\t$seq_name;$seq_length\n";
	$comparison_array[$seqnbr-1] = "$seq_name;$seq_length;"; # start to build up the ordered GMM array
	$seqnbr       ++;
        } #if

    # E.g. SQI total_sequences 66
    if ( $line =~ /^SQI total_sequences (\d+)/ ) {
        $tot_seqs = $1;
        print "\nfound $tot_seqs sequences...\n\n";
        for ( my $i=0; $i<$tot_seqs; $i++ ) {
            $comparison_array[$i] = $seq_names[$i].';'.$seq_length.";";
            }
        } #if
    
    # Example:
    ############ MATRICES ###############
    # typical entry
    # RUN  RANK 0  3  3
    # RUN          BEST_SCORE  8.75
    # RUN INFORMATION_CONTENT  6.44
    # STR    n     S    Lpos   sites       Rpos  seq
    # STR    1  11.40    145   TTTCCTCC    -138  >T07C4.4(spp-1***)
    # STR    2  10.42     42   TTTCCACC     -95  >C53H9.1(rpl-27)
    # STR    3  10.42     32   TTTCCACC    -258  >F32B5.8()
    # ...
    # STR   65   7.72    871   ATTGCTCC    -630  >K11G9.6(mtl-1+)
    # STR   66   6.44    208   TTTGTACC    -374  >T08G5.10(mtl-2+)
    # RUN      ALIGNMENT_CONS  TTTnCTCC
    # RUN        WEIGHTS_CONS  TTTCCTCC
    # CON
    # ALR        1     2     3     4     5     6     7     8
    # ALR   ------------------------------------------------
    # ALR A|     6     0     0     2     2    16     2     7
    # ALR C|     1     5     1    36    57     4    58    50
    # ALR G|     0     0     0     4     0     2     0     3
    # ALR T|    59    61    65    24     7    44     6     6
    # ALR
    # WMR        1     2     3     4     5     6     7     8
    # WMR   ------------------------------------------------
    # WMR A| -0.40 -0.60 -0.75 -0.84 -0.78 -0.11 -0.56 -0.22
    # WMR C| -0.73 -0.30 -0.64  1.51  1.87 -0.55  1.71  1.61
    # WMR G| -0.51 -0.74 -0.42 -0.53 -0.63 -0.62 -0.64 -0.51
    # WMR T|  1.23  1.21  1.40  0.26 -0.06  0.86 -0.10 -0.48
    # WMR    
    #...  
    
    # RUN  RANK 1  5  5
    if ( $line =~ /^RUN\s+RANK\s(\d+)\s+(\d+)\s+(\d+)/ ) {
        $scoring_rank =                 $1;
	$training_run =                 $2;
	$train_iter_best_score_found =  $3;
	$matrix_name =                  "ann_r".$matrix_nbr; # $matrix_name
	$matrix_nbr                     ++;
	print                           $line;
        }
    
    # RUN          BEST_SCORE  8.68
    if ( $line =~     /^^RUN\s+BEST_SCORE\s+(\d+.\d+)/ ) {
        $best_score = $1;
	print         $line;
        }
	
    # RUN INFORMATION_CONTENT  5.70
    if ( $line =~       /^^RUN\s+INFORMATION_CONTENT\s+(\d+.\d+)/ ) {
        $info_content = $1;
	print           "info content: $info_content\n";
        }

    # Site results
    # STR    n     S    Lpos   sites       Rpos  seq
    # STR    1  10.05    576   GAAGAAGT    -925  >Y41D4B.5()
    if ( $line =~ /^STR\s+(\d+)\s+(\S+)\s+(\d+)(\'*)\s+(\w+)\s+(-\d+)\s+>(\S+)/ ) { 
	$n_site_index =            $1;
	$S_wt_mtrx_site_score =    $2;
	$Lpos_Lend_offset =        $3;
	$relend =                  $Lpos_Lend_offset + $width - 1;
	$complement_strand_flag =  $4; # the "'" is only given for a (-)strand hit
	if ( $complement_strand_flag eq "'" ) { $strand = -1; } else { $strand =  1; }
	$sites =                   $5;
	$width =                   length( $sites );
	# the minus sign indicates that Rend offset is in upstream direction
	$Rpos_Rend_offset =        $6;
	$seq_name =                $7; # Y41D4B.5()
        #print                      $line;
	# $hitdata =  "$matrix_name;Lpos_Lend_offset;$relend;$info_content;$strand;";
	$hitdata =                 "$Lpos_Lend_offset;$relend;$info_content;$strand;"; # will add matrix name below
	# to get a matrix name with an alignment consensus attached, we have to hold our breath...
	$hdata{$seq_name} =        $hitdata; # we lose order in the hash, but will get it back by using $seq_name later
    }
    
    # do we want to filter the sorted list of hits and scores to keep just the better scoring hits?
    # ...currently we do this later, by 'thresholding' in the hit plotter
    
    # RUN      ALIGNMENT_CONS  TTTnCTCC
    if ( $line =~ /^^RUN\s+ALIGNMENT_CONS\s+(\S+)/ ) {
        # every time we reach this point, we've just finished another matrix's hits
	$alignment_cons = $1;
	print $line, "\n";

	$matrix_name = $matrix_name."_".$alignment_cons; # add on consensus motif
	
	for ( my $i=0; $i<=$#seq_names; $i++ ) { # go down through the ordered seq list
	    $seq_name =             $seq_names[$i]; # get an ordered seq name
	    # recall a hit data block for this seq -- it will need the matrix name added at its start
	    $this_seq_hitdata =     $hdata{$seq_name}; # this has the incomplete hit data block for this matrix
	    $hdata{$seq_name} =     '';  # empty this hash value so it's ready for the next one
	    $this_seq_hitdata =     $matrix_name.";".$this_seq_hitdata; #now the hit block is complete
	    # so append it to the comparison array
	    $comparison_array[$i] = $comparison_array[$i].$this_seq_hitdata; # append this sequences hit block list
	}
    }
    
    # RUN        WEIGHTS_CONS  TTTCCTCC
    if ( $line =~      /^^RUN\s+WEIGHT_CONS\s+(\S+)/ ) {
        $weight_cons = $1;
	print          $line;	
    }
    
} # READLINE

# Now write out the results array
for ( my $i=0; $i <= $#comparison_array; $i++ ) {
	print        "$comparison_array[$i]\n";
	print OUT    "$comparison_array[$i]\n";
    }

close IN;
close OUT;
exit;