#!/usr/bin/env perl

# transform ChIP-seq eland mappings via hash variable into sorted sga
# usage: eland2sga.pl -f feature -s species eland_file

use strict;
# package to store persistently variables in files [http://search.cpan.org/author/AMS/Storable-2.07/Storable.pm]
use Storable qw (retrieve nstore);
use Getopt::Std;
use File::Spec;

sub Usage {
    die(
"\nusage: eland2sga.pl [-i db-dir] -f feature -a assembly eland_file\n\n"
    );
}

my %opts;

# [TJ] NOTE: The documentation for Getopt::Std (perldoc Getopt::Std) only shows
# examples of colons used with `getopts` (with an 's'), NOT with `getopt`
# (without 's'). In the case of getopt, all options take arguments (as is
# clearly the intention below), but it's not clear that the colon has any effect.
#
getopt( 'asfi', \%opts )
  ; # -a, -s, & -f take arg. Values in %opts, Hash keys will be the switch names

Usage() unless ( $opts{'f'} && ( $opts{'a'} || $opts{'s'} ) );

# hash defining chromosome from SV of chromosomes in current genome assemblies
my $DB;
if ( exists( $opts{'i'} ) ) {
    $DB = $opts{'i'};
}
else {
    $DB = './';
}

# See sga2fps.pl, l. 35
my $nstorage_path = File::Spec->catfile( $DB, "chro_idx.nstorage" );
if ( !( -e $nstorage_path && -r $nstorage_path ) ) {
    print STDERR
"Could not open $nstorage_path: file not found or not readable. Please see option -i.\n";
    Usage();
}
my $chr2AC = retrieve($nstorage_path);

# eventually add older assemblies?

my $species;
if ( exists $opts{'a'} ) {
    $species = $opts{'a'};
}
elsif ( exists $opts{'s'} ) {
    warn "WARNING option -s is deprecated; please use -a.";
    $species = $opts{'s'};
}

my $feature = $opts{'f'};

die("Sorry, $species is not a supported species!\n\n")
  unless ( $chr2AC->{$species} );

my ($score_ref);
while (<>) {
    my $lin = $_;
    chomp $lin;
    my @ar = split( /\t/, $lin );
    $ar[6] =~ s/^.+hromosome\.//;
    $ar[6] =~ s/\.fa$//;

    #    my ($chr)= $ar[6] =~ /hromosome.([0-9RLYX]+).fa/;
    my $chr = 'chr' . $ar[6];
    if ( exists( $chr2AC->{$species}->{$chr} ) ) {

# attribute each tag to its start position and increment counter (multiple column formats)
# >EAS38_3_1_526_564      TGTTCTCTGCATGCTATTTTTTAGTATTGCCGTAA     U0      1       0       0       Mus_musculus.NCBIM37.48.dna.chromosome.15.fa    15876793        R       ..
#.bed:           chr12   119248213       119248236       U0      0       +
        if ( $ar[8] eq 'F' ) {
            $score_ref->{ $chr2AC->{$species}->{$chr} }->{'+'}->{ $ar[7] } += 1;
        }
        elsif ( $ar[8] eq 'R' ) {
            my $pos = $ar[7] + length( $ar[1] );
            $score_ref->{ $chr2AC->{$species}->{$chr} }->{'-'}->{$pos} += 1;
        }
        else {
            warn "parsing error of orientation in following line:\n$lin\n";
        }

# attribute counts to its center position and record in hash (force into numeric)
#        $score_ref->{$chr2AC{$ar[0]}}->{int(($ar[1]+$ar[2])/2)}=eval($ar[3]);
    }
    else {
      #    warn "parsing error of chromosome number in following line:\n$lin\n";
    }
}
#close(FI);

#nstore (\$score_ref, "/scratch/frt/cschmid/EMO-4_bed_ori.nstorage");

my @chr;
my $spec = 'T';    # flag
foreach my $ac ( keys %$score_ref ) {
    unless ( ( $spec eq 'T' ) || ( $spec eq $chr2AC->{'assembly'}->{$ac} ) ) {
        die "The input file ($ARGV) contains entries from both $spec and from ",
          $chr2AC->{'assembly'}->{$ac}, " (chromosome $ac)\n";
    }
    $species = $chr2AC->{'assembly'}->{$ac};
    $chr2AC->{$ac} =~ s/chr//;
    push( @chr, $chr2AC->{$ac} );
}

# order chromosome array as in sga files
my @chrb = sort { $a <=> $b } @chr;
my ( @spl, $sc );
foreach $sc ( 'X', 'Y' ) {
    if ( $chrb[0] eq $sc ) {
        push( @spl, shift(@chrb) );
    }
}
push( @chrb, @spl );

foreach my $chro (@chrb) {

    # foreach position with data, print out counts in sga format:
    #   versioned genome sequence ID
    #   feature type
    #   position
    #   strand (+, -, or 0 for un-oriented features)
    #   tag count (or signal intensity, positive int)

    my @arr =
      keys %{ $score_ref->{ $chr2AC->{$species}->{ 'chr' . $chro } }->{'+'} };
    push @arr,
      keys %{ $score_ref->{ $chr2AC->{$species}->{ 'chr' . $chro } }->{'-'} };

    my %double;    # hash to control for positions with tags in both ori
    foreach my $pos ( sort { $a <=> $b } @arr ) {
        if (
            (
                exists(
                    $score_ref->{ $chr2AC->{$species}->{ 'chr' . $chro } }
                      ->{'+'}->{$pos}
                )
            )
            && ( !$double{$pos} )
          )
        {
            print $chr2AC->{$species}->{ 'chr' . $chro },
              "\t$feature\t$pos\t+\t",
              $score_ref->{ $chr2AC->{$species}->{ 'chr' . $chro } }->{'+'}
              ->{$pos}, "\n";
            $double{$pos} = 'T';
        }
        elsif (
            exists(
                $score_ref->{ $chr2AC->{$species}->{ 'chr' . $chro } }->{'-'}
                  ->{$pos}
            )
          )
        {
            print $chr2AC->{$species}->{ 'chr' . $chro },
              "\t$feature\t$pos\t-\t",
              $score_ref->{ $chr2AC->{$species}->{ 'chr' . $chro } }->{'-'}
              ->{$pos}, "\n";
        }
        else {
            print
"trouble for position $pos on chro 'chr'.$chro: no tags in either orientation??\n";
            die;
        }
    }

    # read length of chromosome from /db/genome/eukaryote.ptr
    my $ac     = $chr2AC->{$species}->{ 'chr' . $chro };
    my $length = $chr2AC->{'length'}->{$ac};
    print "$ac\tEND\t$length\t0\t1\n";

# grep 'Mm/ch' /db/genome/eukaryote.ptr| perl -ane 'if (/^\w+\_\d+\.\d+\s+/){@ar=split/\s+/;print "$ar[0]\tEND\t$ar[2]\t0\t1\n";}'
}

1;

